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Abstract:  The  Training  Range  Environmental  Evaluation  and 
Characterization  System  (TREECS)  is  being  developed  for  the  Army.  This 
system  has  varying  levels  of  capability  to  forecast  the  fate  and  risk  of 
munitions  constituents  (MC)  (such  as  high  explosives  (HE)),  within  and 
transported  from  firing /training  ranges  to  surface  water  and  groundwater. 
The  overall  objective  is  to  provide  Army  environmental  specialists  with  tools 
to  assess  the  potential  for  MC  migration  into  surface  water  and  groundwater 
systems  and  to  assess  range  management  strategies  to  protect  human  and 
environmental  health.  Tier  2  will  consist  of  time-varying  contaminant 
fate/transport  models  for  soil,  vadose  zone,  groundwater,  and  surface  water 
to  forecast  MC  export  from  ranges  and  resulting  concentrations  in  receiving 
waters.  Model  results  can  be  used  to  assess  the  potential  for  surface  water 
and/or  groundwater  MC  concentrations  to  exceed  protective  health 
benchmarks  at  receptor  locations  of  interest. 

The  Tier  2  models  do  not  make  the  highly  conservative  assumptions  of 
steady-state  (time-invariant)  conditions  and  no  MC  loss  or  degradation 
was  used  for  Tier  1.  Thus,  media  concentrations  computed  with  Tier  2 
should  be  closer  to  those  expected  under  actual  conditions.  Media 
concentrations  will  also  reflect  time  phasing  associated  with  time-varying 
MC  loading  conditions  and  transport  arrival  times,  which  can  be  greatly 
extended  for  the  vadose  zone  and  groundwater.  Having  time  as  a 
dimension  in  the  analysis  provides  a  powerful  tool  for  examining  range 
management  strategies  to  promote  attenuation  of  media  concentrations. 

The  information  provided  in  this  report  is  sufficient  to  serve  as  design 
specifications  for  the  development  of  models  and  software  that  will 
comprise  Tier  2  of  TREECS.  The  details  of  the  Tier  2  soil  model 
formulations  provided  herein  can  also  help  serve  as  documentation  for 
that  model.  All  components  will  be  packaged  within  a  user-friendly  PC 
client -based  application  with  an  emphasis  on  ease-of-use. 


DISCLAIMER:  The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products. 
All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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meters 
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flow,  m 

Li,k  MC  residue  mass  loading  for  constituent  i  for  year  k,  g/yr 

LOj,k  percent  of  low  order  detonations  for  munitions  item  /  for 

year  k 

L(t)  time-varying  solid  phase  MC  mass  loading,  g/yr 

/  length  of  a  right  cylindrical  MC  particle,  cm 

Ma  MC  vapor  mass  in  air  space  within  soil  matrix,  g 

Md  MC  mass  dissolved  in  soil  matrix  pore  water,  g 

Mdis(t)  particle  mass  loss  over  time  due  to  dissolution,  g 

My  mass  of  constituent  i  in  munitions  item  /  delivered  to  impact 

area,  g/item 

Mns  non-solid  phase  MC  mass,  g 

Mp  MC  mass  adsorbed  to  soil  particles,  g 
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Ms  solid  phase  MC  mass,  g 

Msoii  mass  of  soil  in  the  AOI,  g 

Mi  molecular  weight  of  the  MC  of  interest 

M2  molecular  weight  of  carbon  disulfide,  76  g/mole 

m  total  solid  phase  MC  mass  for  all  solid  phase  MC  particles,  g 

m  MEPAS  groundwater  model  coefficient  dependent  on  soil 

type  where  m  =  2b  +  3,  dimensionless 
rru  MC  solid  phase  mass  for  particle  i,  g 

N  average  number  of  rainfall  events  per  year,  events/yr 

Nj,k  number  of  munitions  item  /  fired  for  year  k 

n  total  number  of  munitions  items  used  at  AOI 

P  average  annual  precipitation,  m/yr 

Pr  average  annual  rainfall,  m/yr 

Pt  average  annual  total  precipitation,  cm/yr  or  m/yr  depending 

on  use 

Q  average  annual  water  flow  rate  through  the  water  body  for 

RECOVERY  model,  m3/yr 

Q  constant  background  stream  flow  rate  (e.g.,  annual  mean 

flow)  at  the  head  of  the  reach  (without  any  flow  from  the 
AOI)  for  CMS,  ms/yr 

Qw  water  flow  rate  due  to  net  infiltration,  or  percolation,  from 

soil  (groundwater  recharge  into  vadose  zone),  m3/yr 
qbv  volume  of  bed  material  load  per  unit  width  in  Duboy’s 

equation,  ft2/sec 

qbv *  dimensionless  volumetric  unit  sediment  discharge 

qr  average  annual  surface  water  runoff  rate,  m/yr 

qw  average  annual  Darcy  water  infiltration  rate,  m/yr 

R  retardation  factor,  dimensionless 

RedQ  event -based  runoff  mass  removal  rate  of  pore  water  due  to 

rain-induced  ejection,  g/sec 

Ru  universal  gas  constant  =  8.206  E-5  atm0m3/g-mol  °K 

R/s  retardation  factor  in  the  saturated  zone  for  the  MEPAS 

groundwater  model,  dimensionless 

Rfu  retardation  factor  for  unsaturated  zone  or  vadose  zone  for 

the  MEPAS  groundwater  model,  dimensionless 
S  TSS  concentration  for  CMS,  g/ms 

S  water  surface  slope,  ft/ft 

Sw  TSS  input  variable  name  for  RECOVERY,  g/m°- 

SYMj,k  percent  of  sympathetic  detonation  of  duds  for  munitions 
itemj  for  year  k 
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T  time  averaging  interval  for  rainfall  event,  or  the  rainfall  event 

duration,  sec 

T1/2  half-life  due  to  decay  or  degradation,  yr 

T  average  or  typical  rainfall  event  duration,  time 

T  water  and  sediment  mean  temperature  for  CMS,  °C 

Ta  ambient  absolute  temperature  of  the  soil  matrix,  °K 

Tp  total  time  period  of  simulation  for  CMS,  yr 

Tsoii  average  annual  temperature  of  the  soil  matrix,  °C 

Tw  ambient  water  temperature  for  dissolution,  °C 

t  time,  yr  and  sec 

td  rain  drop  arrival  interval  for  explosives  dissolution,  sec 

tgo  time  required  to  reached  90  percent  mass  loss  from 

dissolution,  yr 

U  stream  mean  velocity  for  CMS,  m/day  and  m/sec 

Uw  mean  wind  speed  for  RECOVERY  model,  m/sec 

u  pore  water  velocity  in  horizontal  (X)  direction  for  the 

MEPAS  groundwater  model,  cm/sec 
uc  critical  flow  velocity  for  particle  movement,  ft/sec 

u*  shear  velocity  of  the  flow  for  CMS,  m/sec 

V  surface  soil  compartment  volume  of  AOI,  m3 

Va  volume  of  air  in  the  AOI  soil,  m3 

Vb  active  sediment  layer  burial  rate  for  CMS  and  REOVERY 

models,  m/day  and  m/sec 

Vd  mass  transfer  rate  across  the  sediment -water  interface 

resulting  from  diffusion  of  the  dissolved  constituent  for  CMS 
and  RECOVERY  models,  m/day  and  m/sec 
Vd  Darcy  velocity  of  groundwater  flow,  L/T 

Vi  water  layer  volume  coating  a  spherical  explosive  particle, 

cm3 

Vr  sediment  resuspension  rate  for  CMS  and  RECOVERY 

models,  m/day  and  m/sec 

Vs  suspended  solids  settling  rate  for  CMS  and  RECOVERY 

models,  m/day  and  m/sec 
Vw  volume  of  water  in  the  AOI  soil,  m3 

W  additional  constant  external  loading  rate  of  constituent  for 

RECOVERY  model,  kg/yr 

W  mean  wind  speed  at  10  m  above  surface  for  CMS,  m/sec 

w  pore  water  velocity  in  the  vertical  (Z)  direction  of  the  vadose 

zone  for  the  MEPAS  groundwater  model,  cm  S’1 
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Wf  width  of  the  AOI,  or  dimension  perpendicular  to  the 

groundwater  flow,  m 

X  longitudinal  coordinate  for  the  MEPAS  groundwater  model 

X  Usage  location,  i.e.,  distance  downstream  from  the  upstream 

boundary  to  the  location  of  interest  for  CMS,  km 
x  longitudinal  distance  from  the  center  of  the  AOI  to  the 

receptor  well  for  the  MEPAS  groundwater  model,  cm 
x  downstream  distance  along  stream  for  CMS,  m 

x  particle  shape  index  for  solid  phase  MC  particles  that  is  input 

for  the  soil  model;  l  for  cylindrical  and  2  for  spherical, 
dimensionless 

Y  lateral  coordinate  for  the  MEPAS  groundwater  model 

y  lateral  distance  from  the  groundwater  plume  centerline  to 

the  receptor  well  for  the  MEPAS  groundwater  model,  cm 
YHOj,k  percent  yield  of  munitions  item  j  due  to  high  order 

detonation  for  year  k 

Ywj,k  percent  yield  of  munitions  item  j  due  to  low  order  detonation 

for  year  k 

YsYMj,k  percent  yield  of  munitions  item  j  due  to  sympathetic 

detonation  for  year  k 

Z  vertical  coordinate  for  the  MEPAS  groundwater  model 

Za  aquifer  thickness,  L 

Zb  surficial  soil  layer  thickness  of  AOI,  m 

Zv  vadose  zone  thickness,  L 

z  distance  along  the  vertical  coordinate  or  depth  of  the 

receptor  well  below  the  water  table  for  the  MEPAS 
groundwater  model,  cm 

z  depth  of  mixed  sediment  layer  for  RECOVERY  model,  m 

a  average  specific  surface  area  of  the  solid  phase  MC  mass, 

cm2/g 

ax  longitudinal  dispersivity  for  the  MEPAS  groundwater  model, 

cm 

ay  lateral  dispersivity  for  the  MEPAS  groundwater  model,  cm 

az  vertical  dispersivity  for  the  MEPAS  groundwater  model,  cm 

a0  initial  particle  radius  for  dissolution,  cm 

P  consolidating  parameter  used  in  calculating  rainfall- 

extraction  (see  Equation  27),  sec1 
pm  MC  mass  dissolution  flux  rate,  g/cm2-sec 
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At  time  step  in  the  solution  procedure  for  the  Tier  2  soil  model, 

yr 

C,  MC  dissolution  mass  transfer  rate  or  solid  phase  surface 
removal  rate,  cm/sec 

y  specific  weight  of  water,  62.4  lb/ft3 

y  consolidating  parameter  for  dissolution  parameters, 

r  =  P,aCs 

ys  specific  weight  of  sediment,  dimensionless 

k  a  computed  parameter  used  in  rainfall-extraction 

calculations  (see  Equation  33),  dimensionless 
A  degradation  or  decay  rate  for  the  MEPAS  groundwater 

model,  sec-1 

Aa  degradation  rate  for  the  aqueous  adsorbed  MC,  yr1 

Ai  degradation  rate  for  the  liquid  (aqueous  dissolved)  MC,  yr1 

<f>  soil  porosity  or  ratio  of  void  volume  to  total  volume;  water 

content  when  water-saturated  and  porosity  of  the  sediment 
bed  for  CMS,  fraction 

(f)(,  effective  porosity  for  the  MEPAS  groundwater  model, 

fraction 

(jhn  mixed  sediment  layer  porosity  for  RECOVERY  model, 

fraction 

<jh  deep  sediment  porosity  for  RECOVERY  model,  fraction 

pb  soil  dry  bulk  density,  g/ml  or  kg/L 

ps  dry  sediment  particle  density,  kg/L 

psm  solid  phase  constituent  mass  density,  g/cms 

6  an  empirical  coefficient  to  account  for  temperature  effects  on 

dissolution  rate  due  to  changes  in  solubility 
Of  soil  field  capacity,  fraction 

6w  soil  volumetric  moisture  content  or  ratio  of  water  volume  to 

total  volume,  fraction 

zb  critical  shear  stress  for  particle,  lb/ft2 

To  shear  stress  of  the  water  flow,  lb/ft2 

r w  water  residence  time  for  RECOVERY  model,  yr 

r*  dimensionless  Shields  parameter 

coo  Rubey’s  clear-water  sediment  fall  velocity,  L/T 
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1  Introduction 

Background 

The  Training  Range  Environmental  Evaluation  and  Characterization 
System  (TREECS)  is  being  developed  for  the  Army.  The  system  has  varying 
levels  of  capability  to  forecast  the  fate  of  munitions  constituents  (MC)  (such 
as  high  explosives  (HE))  within  and  transported  from  firing/training  ranges 
to  surface  water  and  groundwater.  The  overall  objective  is  to  provide  Army 
environmental  specialists  with  tools  to  assess  the  potential  for  MC 
migration  into  surface  water  and  groundwater  systems  and  to  assess  range 
management  strategies  to  protect  human  and  environmental  health. 

TREECS  will  be  accessible  from  the  World  Wide  Web  and  will  initially  have 
two  tiers  for  assessments.  Tier  l  (Dortch  et  al.  2009)  consists  of  screening- 
level  methods  that  require  minimal  data  input  requirements  and  can  be 
easily  and  quickly  applied  by  Army  environmental  staff  to  assess  the 
potential  for  MC  migration  to  surface  water  and/or  groundwater  MC  with 
concentrations  exceeding  protective  health  benchmarks.  Various 
simplifying  assumptions,  such  as  steady-state  conditions  and  no  system 
losses,  such  as  degradation,  are  made  to  provide  conservative  or  worst-  case 
estimates  with  Tier  1.  If  a  potential  concern  is  indicated  by  a  Tier  1  analysis, 
then  there  would  be  cause  to  proceed  to  Tier  2  to  obtain  a  more  definitive 
assessment. 

Tier  2  assessment  methods  will  require  more  detailed  site  data,  and  will 
require  more  knowledge  and  skill  to  apply.  However,  these  methods  can  be 
applied  by  environmental  staff  who  have  a  cursory  understanding  of  multi- 
media  fate  and  transport.  The  Tier  2  approach  will  allow  time-varying 
analyses  and  system  losses.  Such  analyses  should  provide  more  accurate 
predictions  with  generally  lower  concentrations  due  to  mediating  effects  of 
transport  phasing/dampening,  natural  attenuation,  and  degradation.  Tier  2 
should  prove  even  more  valuable  than  Tier  1  for  range  management 
assessments  given  the  fact  that  it  includes  the  time  domain;  thus  allowing 
assessment  of  the  effectiveness  of  range  use  rotation  over  time  as  a 
management  option.  Both  Tiers  1  and  2  focus  on  contaminant  stressors  and 
human  and  ecological  health  end  point  metrics. 
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Objective 

The  objective  of  this  report  is  to  describe  the  approach  and  formulations 
that  will  be  used  for  Tier  2  of  TREECS.  The  information  provided  in  this 
report  is  sufficient  to  serve  as  design  specifications  for  the  development  of 
models  and  software  that  will  comprise  Tier  2.  The  details  of  the  soil 
model  formulations  provided  herein  can  also  serve  as  documentation  for 
that  model. 

Scope 

The  Army’s  Operational  Range  Assessment  Program  (ORAP)  is  being 
conducted  in  two  phases:  a  Qualitative  Assessment  (Phase  I)  and  a 
Quantitative  Assessment  (Phase  II).  The  results  of  Phase  I  placed  ranges 
into  one  of  two  categories:  unlikely  and  inconclusive.  A  range  categorized 
as  “unlikely”  requires  no  further  action,  and  it  is  placed  into  a  five-year 
review  cycle.  Ranges  categorized  as  “inconclusive”  will  require  a  follow-up 
Phase  II  assessment. 

The  Army  completed  all  ORAP  Phase  I  assessments  at  the  end  of  FY  09. 
Phase  II  assessments  commenced  in  FY  10  after  the  completion  of  all 
Phase  I  assessments.  Given  this  timetable  and  the  more  comprehensive 
needs  and  challenges  of  the  Phase  II  assessments,  ERDC,  in  consultation 
with  the  Army  Environmental  Command  (AEC),  concluded  that  the  Tier  1 
and  Tier  2  tools  of  TREECS  could  help  support  needs  within  the  Phase  II 
ORAP.  Additionally,  these  tools  should  help  address  range  management 
issues,  such  as  how  to  operate  a  range  to  avoid  problems  with  off-site 
migration  of  MC. 

In  general,  the  Tier  1  and  2  tools  of  TREECS  should  provide  the  following 
information  to  address  the  needs  of  the  Army: 

1.  Given  range  use,  estimate  the  MC  residue  mass  loading  rate  to  the  range  as 
mass  (M)  per  time  (T),  M/T. 

2.  Given  the  mass  loading  rate  of  MC,  estimate  the  soil  concentration  on  the 
range  area  of  interest  and  the  mass  fluxes  (M/T)  off  the  range  to  other 
media  (e.g.,  surface  water  and  groundwater). 

3.  Given  the  mass  fluxes  to  other  media,  estimate  the  media  concentrations 
at  points  of  interest  off-range. 

4.  Given  the  media  concentrations  at  points  of  interest  off-range,  determine  if 
protective  health  benchmark  concentrations  are  exceeded. 
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This  information  can  help  environmental  specialists  determine  whether 
there  is  potential  for  a  particular  range  to  develop  into  a  human  health  or 
environmental  hazard.  Furthermore,  range  use  strategies  could  be 
adjusted  or  managed  in  an  effort  to  promote  range  sustainment  while 
satisfying  protective  health  requirements.  Having  an  estimate  of 
concentration  versus  time  in  groundwater  and  surface  water  will  provide 
additional  information  that  could  be  used  in  the  design  of  sampling  and 
monitoring  strategies. 

Both  Tiers  l  and  2  can  address  the  above  needs,  but  the  approach  to  each 
is  different.  The  primary  difference  between  the  two  is  that  Tier  1  is  time- 
invariant  or  steady-state,  and  Tier  2  is  time-varying.  Including  time- 
varying  conditions  in  the  analysis  substantially  increases  the  level  of 
accuracy  in  the  forecasts  at  the  cost  of  needing  more  data  and  knowledge. 
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2  Approach 

The  approach  taken  for  Tier  2  is  similar  to  Tier  1  in  many  aspects, 
especially  with  regard  to  spatial  dimensionality.  However,  there  are  far 
fewer  simplifying  assumptions  made  for  Tier  2,  as  compared  with  Tier  1, 
since  time  is  included  in  the  analysis  approach.  The  conceptual  site  model 
and  the  assumptions  that  will  be  made  for  Tier  2  are  summarized  in  this 
chapter.  Descriptions  of  the  models  are  provided  later  in  the  chapters 
describing  each  model. 

Conceptual  Model  Description 

For  land-based  firing  ranges,  four  media  are  considered  when  determining 
MC  fate:  soil,  vadose  zone,  groundwater  or  aquifer,  and  surface  water 
(including  surface  water  sediments).  Potential  air  concentrations  of  MC 
are  a  limited,  short-duration,  local  issue  and  are  not  considered  for  range 
sustainment.  A  conceptual  site  model  (CSM)  for  range-generated  MC  is 
shown  in  Figure  1.  MC  residue  loading  first  enters  the  range  soil.  MC  can 
move  from  soil  to  surface  water  via  rainfall-extracted  runoff  and  erosion 
and  from  soil  to  the  vadose  zone  via  infiltration  or  leaching.  MC  can  then 
percolate  through  the  vadose  zone  into  a  receiving  aquifer.  Aquifer  and 
surface  water  concentrations  of  MC  depend  on  the  location  of  interest 
within  the  receiving  water,  relative  to  the  point  of  MC  influx.  Receiving 
water  concentrations  can  be  computed  and  compared  with  benchmarks  for 
compliance.  There  can  also  be  pathways  from  groundwater  to  surface 
water  and  vice  versa;  however,  only  the  pathway  from  groundwater  to 
surface  water  is  considered  important  in  Tier  2.  There  can  be  a  pathway 
for  interflow  from  soil  to  surface  water  in  some  situations  as  well. 

Interflow  is  horizontal  water  movement,  rather  than  vertical  infiltration, 
within  the  soil  or  vadose  zone.  This  feature  will  be  included  in  Tier  2  and 
will  be  discussed  in  Chapter  8. 

The  main  differences  between  the  Tier  1  CSM  (Dortch  et  al.  2009)  and  the 
Tier  2  CSM  shown  herein  is  that  Tier  2  includes  the  vadose  zone  and  a  link 
from  groundwater  to  surface  water.  The  link  from  groundwater  to  surface 
water  can  also  be  added  later  to  Tier  1.  Also,  Tier  2  includes  system  losses, 
including  degradation  and  volatilization,  and  solid  phase  mass  dissolving 
to  non-solid  phase  mass,  whereas  Tier  1  does  not. 
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Figure  1.  Tier  2  conceptual  model  schematic. 


With  the  inclusion  of  time-varying  conditions  in  Tier  2,  it  is  necessary  to 
include  the  vadose  zone  since  mass  transit  time  through  the  vadose  zone 
can  be  quite  long  for  some  systems.  There  can  also  be  mass  loss  within  the 
vadose  zone  due  to  degradation. 

The  soil  model  will  compute  time-varying  soil  concentrations  and  mass 
export  fluxes  (M/T)  for  erosion,  rainfall-extracted  runoff,  and  infiltration. 
A  portion  of  the  infiltration  flux  can  be  designated  as  interflow  with  export 
to  surface  water  rather  than  the  vadose  zone. 

The  vadose  zone  model  will  use  the  net  infiltration,  or  leached  mass  influx 
rate,  to  compute  the  time-varying  mass  flux  (M/T)  entering  the  aquifer. 
The  aquifer  model  will  use  the  vadose  zone  mass  flux  and  the  receptor 
location  to  compute  the  time-varying  groundwater  concentration  at  the 
receptor  location.  Sorption  and  degradation  can  exist  within  the  aquifer. 

Mass  flux  imported  by  surface  water  can  include  soil  compartment  exports 
of  runoff  and  erosion,  plus  interflow  and  groundwater  discharge  to  surface 
water.  The  surface  water  models  will  use  the  time-varying  mass  influx 
imports  to  compute  the  time-varying  surface  water  and  sediment 
concentrations.  There  can  be  sorption,  degradation,  volatilization,  mass 
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transfer  between  the  water  column  and  sediment  bed,  and  sediment  burial 
in  the  bed.  The  location  of  the  receptor  will  not  be  required  for  one  of  the 
surface  water  models,  but  it  is  required  for  another  surface  water  model 
that  will  be  included  for  optional  use.  The  latter  surface  water  model 
computes  concentration  at  the  receptor  location  as  described  later.  The 
models  for  each  medium  and  their  inputs  and  output  are  described  later  in 
this  report. 

The  proof-of-concept  testing  of  Tier  l  (Dortch  et  al.  2010)  revealed  that 
certain  conditions  needed  to  be  imposed  on  the  CSM  and  its  implemen¬ 
tation  within  TREECS.  These  conditions  apply  to  Tier  1  and  2,  and  they 
include  the  following:  each  application  will  be  for  a  single  area  of  interest 
(AOI),  and  each  AOI  constitutes  a  separate  application;  each  AOI  will  have 
no  more  than  one  aquifer  and  one  surface  water-body  receiving  loadings 
from  the  AOI  being  assessed. 

When  an  impact  area  or  AOI  is  within  a  single  drainage  basin,  or  sub-basin, 
of  a  watershed,  designation  of  the  AOI  is  relatively  straightforward;  the 
impact  area  should  constitute  a  single  AOI  application.  However,  when  an 
impact  area  is  rather  large  and  overlaps  more  than  one  sub-basin,  it  may  be 
more  accurate  to  split  the  impact  area  into  multiple  AOI.  Such  an  example  is 
shown  in  Figure  2  where  the  impact  area  overlaps  two  drainage  basins.  This 
example  analysis  could  be  conducted  for  two  AOI,  A  and  B.  AOI  A  would  use 
Lake  A  for  the  receiving  water  analysis,  and  AOI  B  would  use  Lake  B.  The 
difficulty  with  conducting  an  analysis  like  this  is  distributing  the  MC 
loadings  between  AOI  A  and  B.  Many  firing  ranges  often  fire  into  the  same 
impact  area.  Range  records  usually  denote  the  types  and  numbers  of  muni¬ 
tions  fired  from  each  range,  but  they  do  not  describe  where  the  projectiles 
hit.  Thus,  it  is  not  possible  to  assign  the  AOI  receiving  the  projectiles  and 
their  residue.  One  approach  could  be  to  distribute  the  residue  loadings 
according  to  the  fraction  of  total  impact  area  that  each  AOI  constitutes. 

Basic  Assumptions 

As  with  Tier  1,  each  AOI  will  be  treated  as  a  single  fully  mixed  compartment; 
thus,  soil  concentrations  are  assumed  to  be  uniform  throughout  the  upper 
soil  horizon  within  the  AOI.  As  an  example,  the  primary  impact  zone  of  a  set 
of  ranges  will  be  treated  as  a  single  homogenous  area.  Although  an  AOI  is, 
in  fact,  heterogeneous,  treating  a  heterogeneous  AOI  as  homogeneous  is  not 
a  fatal  assumption  because  the  total  MC  source  mass  loading  to  soil  is  the 
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Figure  2.  Example  of  an  impact  area  that  overlaps  two  drainage  basins. 


driving  variable  for  export  flux, not  MC  concentration  in  soil.  Although  the 
soil  concentration  of  MC  depends  on  heterogeneities  and  even  the  size  of 
the  source  area  for  a  homogeneous  site,  the  fluxes  or  export  rate  of  mass 
from  the  source  area  to  other  media  does  not  depend  on  the  area  of  the 
source  zone  or  the  source  concentration.  Of  course,  a  large  impact  area  with 
source  mass  clusters  that  differ  substantially  could  be  treated  as  multiple 
sources  with  each  assessed  individually.  Predicted  receiving  media 
concentrations  will  be  more  conservative  (i.e.,  higher  concentrations) 
nonetheless  if  multiple  source  areas  and  the  respective  masses  are 
combined. 

The  loading  of  MC  into  the  AOI  will  be  estimated  from  the  numbers  and 
types  of  munitions  used  on  an  annual  basis.  Additionally,  the  percentages 
for  low  order  detonations,  duds,  and  sympathetic  detonations  will  be 
assigned  for  each  type  of  munitions.  Sympathetic  detonations  are  duds  that 
undergo  a  low  order  detonation  due  to  a  nearby  high  order  detonation.  The 
munitions  yield,  or  percent  of  exploded  energetic  mass,  must  also  be 
assumed  and  assigned  for  each  munitions  type  and  its  associated  detonation 
type.  Initially,  Tiers  l  and  2  of  TREECS  will  address  only  impact  areas  and 
will  not  include  firing  points,  but  firing  points  could  potentially  be  added 
later.  Additionally,  the  user  will  have  the  option  to  specify  a  known  loading 
rate  of  MC  that  may  have  been  estimated. 
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Time-varying  MC  loadings  and  receiving  media  response  can  be  assessed 
with  Tier  2  since  the  steady-state  assumption  used  for  Tier  1  will  not  be 
made.  Allowing  for  time-varying  conditions  increases  model  input 
requirements,  but  it  also  greatly  increases  the  amount  of  output  information 
available  and  can  affect  how  receiving  waters  respond  to  range  loadings  and 
export.  For  example,  if  a  loading  of  MC  occurs  for  a  ten -year  period  and 
then  ceases,  it  may  take  another  10  years  after  the  loading  ceases  for 
groundwater  concentration  to  peak,  and  the  arrival  time  and  concentration 
level  of  the  peak  varies  with  the  distance  of  the  well  from  the  source  area 
and  environmental  factors.  These  variations  can  attenuate  the  groundwater 
concentrations  to  levels  that  can  be  substantially  lower  than  would  occur  for 
continuous,  long-term,  steady-state  MC  loadings.  Consequently,  it  is 
expected  that  receiving  water  concentrations  will  be  generally  lower  for  a 
Tier  2  analysis  compared  with  a  Tier  1  analysis. 

The  input  requirements  will  be  greater  for  Tier  2  than  for  Tier  1.  System 
loss  processes,  such  as  degradation  and  volatilization,  require  input 
parameters.  With  a  time-varying  soil  model,  it  will  be  possible  to  include 
corrosion/dissolution  rates  that  drive  the  transformation  of  a  solid  phase 
MC  mass  to  the  more  mobile  aqueous  phase.  There  are  input  parameters 
that  will  be  required  for  dissolution  processes.  In  addition,  it  will  be 
necessary  to  input  the  initial  soil  concentration  of  MC  when  applying 
Tier  2.  Tier  2  input  requirements  are  discussed  in  greater  detail  later  in 
this  report. 

The  receiving  surface  water  can  be  a  stream,  pond,  lake,  wetland,  or  any 
type  of  surface  water.  It  will  be  assumed  that  any  erosion  mass  trapped 
within  the  AOI  drainage  area  prior  to  export  from  the  AOI  can  described 
with  a  sediment  delivery  ratio,  or  SDR,  as  discussed  in  Chapter  4.  The  use 
of  the  SDR  is  optional  and  will  be  calculated  within  the  soil  model  and 
applied  to  eroded  export  fluxes  computed  by  the  soil  model.  It  is  assumed 
that  there  are  no  conveyance  losses  or  trapping  of  MC  mass  between  the 
AOI  and  receiving  surface  water-body  for  AOI  exports  to  surface  water; 
this  provides  a  conservative,  worst-case  scenario  while  greatly  reducing 
model  complexity  and  data  input  requirements. 

Transport  through  the  vadose  zone  from  soil  to  the  groundwater  will  be 
treated  as  a  one-dimensional  (lD),  vertical  transport  process.  Lateral 
transport  is  not  considered,  and  vertical  transport  is  assumed  to  be 
uniform  horizontally  over  the  AOI.  The  aquifer  transport  is  described  as 
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unidirectional  (lD),  horizontal  flow  (advection)  with  three-dimensional 
(3D)  dispersion  along  the  flow  axis.  Degradation  in  the  vadose  zone  and 
aquifer  is  assumed  to  be  first  order.  Reversible,  linear,  equilibrium 
sorption  is  assumed  in  both  media. 

Two  options  will  be  provided  for  modeling  MC  fate  in  surface  water.  One 
option  is  the  same  model  used  in  Tier  1  for  surface  water,  which  is  the 
RECOVERY  model.  The  second  option  is  the  Contaminant  Model  for 
Streams  (CMS).  The  RECOVERY  model  treats  the  surface  water  column  as 
a  single,  fully  mixed  compartment,  which  is  zero-dimensional  (oD).  The 
water  column  is  underlain  with  a  lD  series  of  vertical  layers  with  varying 
properties  and  concentrations.  The  CMS  treats  the  surface  water  as  a  series 
of  nodes  where  concentrations  are  uniform  at  a  node  (i.e.,  uniform  over  the 
depth  and  width  of  the  water)  but  vary  from  node  to  node  or  along  the 
length  of  the  water  body  (i.e.,  in  the  longitudinal  direction).  Each  node 
includes  a  benthic  sediment  compartment  in  addition  to  the  water  column, 
and  there  can  be  mass  transferred  between  water  and  sediment;  but  no 
mass  transfer  is  allowed  between  sediment  compartments  in  adjacent  nodes 
(i.e.,  oD  at  a  node).  Sediment  concentrations  can  also  vary  from  node  to 
node  in  the  longitudinal  direction.  The  RECOVERY  model  is  best  suited  for 
standing  or  pooled  water,  such  as  pond  and  lakes.  CMS  is  best  suited  for 
freely  moving  water,  such  as  streams  and  rivers.  Either  model  can  be 
applied  to  any  type  of  surface  water  body,  but  the  user  should  understand 
the  ramifications  of  the  basic  assumptions:  i.e.,  oD  water  column  and  lD 
sediment  representation  in  RECOVERY  and  lD  water  column  in  CMS.  The 
CMS  is  best  suited  for  tidal  streams  as  discussed  in  Chapter  6. 

Like  Tier  1,  Tier  2  will  use  steady-state  hydrologic  inputs  to  compute  model 
fate  processes.  The  long-term  average  rates  for  annual  soil  erosion,  rainfall 
and  precipitation,  and  infiltration  will  be  used  for  the  soil  model.  Similarly, 
input  rate  parameters  for  fate  processes  are  assumed  to  be  constant  over  the 
simulation.  The  use  of  average  annual  hydrology  and  rate  parameters  does 
not  preclude  employing  a  time-varying  contaminant  fate  model  since  the 
contaminant  mass  balance  equations  are  time-varying.  Time-varying 
concentrations  and  export  fluxes  will  result  even  with  constant  input 
parameter.  However,  the  time-varying  results  will  approach  steady-state  if 
the  MC  loadings  are  constant  in  time.  The  use  of  constant  input  parameter 
greatly  reduces  model  complexity  and  increases  model  ease-of-use. 
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The  formulations  or  descriptions  for  the  models  used  within  Tier  2  of 
TREECS  are  described  in  Chapters  3  through  6.  These  models  include:  MC 
residue  mass  loading  to  soil,  MC  fate  and  export  in  soil,  MC  fate/transport 
in  vadose  zone,  MC  fate/transport  in  aquifer,  and  MC  fate/transport  in 
surface  water  and  surface  water  sediment.  Other  modeling  assumptions 
are  addressed  within  each  of  these  chapters. 
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3  MC  Residue  Mass  Loading  Model 

The  MC  residue  mass  loading  model  for  Tier  2  will  be  similar  to  the  one  for 
Tier  1;  the  primary  difference  is  that  munitions  use  can  be  varied  on  an 
annual  basis  rather  than  being  constant  over  time.  For  each  munitions  item 
used,  the  user  first  selects  the  munitions  identification  using  the  munitions 
type  and  the  Department  of  Defense  identification  code  (DODIC)  or 
National  Stock  Number  (NSN).  For  each  item,  the  user  then  provides  the 
following  for  each  year  of  input: 

•  The  year  that  the  item  was  used 

•  The  number  fired  per  year 

•  The  percent  of  duds 

•  The  percent  of  low  order  detonations 

•  The  percent  yield  (portion  of  MC  used  up  when  munitions  explode)  for 
low  order  detonations 

•  The  percent  of  duds  that  are  sympathetically  exploded  by  another 
detonation 

•  The  percent  yield  for  sympathetic  detonations 

•  The  percent  yield  for  high  order  detonations. 

TREECS  will  include  help  files  for  setting  the  above  percentages.  The 
percent  of  high  order  detonations  is  calculated  from  too  percent  minus  the 
percentages  of  duds  and  low  orders,  and  the  calculated  value  is  entered  in 
the  input  table.  The  user  cannot  change  this  displayed  value  except  for 
changing  either  the  percentage  of  duds  or  low  orders. 

The  item  usage  per  year  is  stepped;  i.e.,  the  input  numbers  remain  constant 
until  the  next  update  year  in  the  input  table.  An  example  of  the  input  screen 
for  the  Tier  2  MC  residue  mass  loading  model  is  shown  in  Figure  3.  In  this 
example,  the  usage  for  DODIC  C445  is  500  firings  per  year  from  1950 
through  1954.  In  1955,  the  firing  of  this  item  increased  to  1000  per  year. 
Notice  that  the  other  input  parameters  can  also  vary  annually.  If  only  one 
line  of  annual  input  is  supplied  within  the  munitions  usage  table  for  each 
munitions  item,  then  the  residue  mass  loading  would  be  constant  over  time, 
or  time-invariant,  as  in  Tier  1  TREECS. 
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Figure  3.  Example  of  input  screen  for  Tier  2  MC  residue  mass  loading  model. 

In  order  to  conduct  the  estimates  of  MC  residue  loading,  it  is  necessary  to 
first  know  the  amount  of  MC  mass  in  each  munitions  item  that  is  delivered 
to  the  impact  area.  This  information  can  be  obtained  from  the  Munitions 
Items  Disposition  Action  System  (MIDAS)  (https://midas.dac.armv.mil/)  based  on 
DODIC  or  NSN.  However,  extraction  of  information  from  MIDAS  can  be 
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slow  and  tedious.  A  utility  has  been  developed  for  automatically  pulling 
this  information  into  the  TREECS  application  using  the  DODIC  or  NSN. 
This  utility  requires  a  special,  processed  subset  database  of  MIDAS  that 
was  developed  by  the  Defense  Ammunition  Center,  DAC.  The  MIDAS 
subset  database  presently  includes  data  for  164  DODIC  items,  but  this 
database  can  and  should  be  expanded  to  include  more  items. 


Once  the  MC  mass  delivered  to  the  impact  area  is  known  for  each 
munitions  item  used,  the  calculation  of  residue  mass  loadings  is  a 
straightforward  summation.  The  MC  residue  mass  loading  for  constituent 
i  for  year  k,  U,k,  (g/yr)  is  computed  as  follows, 


L 


i,k 


where 

DUDj,k 

HOj,k 


LOj,k 
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Yhojm 
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LO  100  -Y 


)  +  HO  (lOO  Y  )  +  DUD  SYM  (lOO 
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(1) 


=  percent  of  duds  for  munitions  item  j  for  year  k 
=  percent  of  high  order  detonations  for  munitions  itemj  for  year 
k 

=  percent  of  low  order  detonations  for  munitions  itemj  for  year 
k 

=  mass  of  constituent  i  in  munitions  item  /  delivered  to  impact 
area,  g/item 

=  number  of  munitions  item  j  fired  for  year  k 
=  total  number  of  munitions  items  used  at  AOI 
=  percent  of  sympathetic  detonation  of  duds  for  munitions  item  / 
for  year  k 

=  percent  yield  of  munitions  item  j  due  to  high  order  detonation 
for  year  k 

=  percent  yield  of  munitions  item  j  due  to  low  order  detonation 
for  year  k 

=  percent  yield  of  munitions  item  j  due  to  sympathetic 
detonation  for  year  k 


Thus,  the  number  of  items  fired,  the  detonation  type  percentages,  and  the 
percent  yields  can  vary  on  an  annual  basis.  As  stated  previously,  the  percent 
of  low  order  detonations,  duds,  and  high  order  detonations  must  sum  to 
too  percent  for  each  year  of  input.  Additionally,  annual  input  values  must 
be  provided  for  each  munitions  item  designated  as  used  at  the  AOI.  There  is 
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fairly  good  guidance  on  dud  and  low  order  rates  (Concurrent  Technologies 
Corporation  2004)  for  detonations,  but  there  is  only  limited  information  on 
percent  yield  (Gerald  et  al.  2007).  There  is  currently  no  guidance  on 
percentage  of  sympathetic  detonations.  The  residue  mass  loading  model  will 
also  include  the  capability  to  specify  a  user-defined  or  known  mass  loading 
rate  for  each  constituent  rather  than  using  the  rates  computed  based  on 
munitions  usage. 
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4  Time-Varying  Soil  Model 

Consider  a  compartment  of  surface  soil  of  thickness  or  depth  Zb  (m)  over  a 
given  area  of  interest  (AOI)  with  surface  area  A  (m2).  Assume  that  the  MC 
concentration  varies  with  time  but  is  homogeneous  or  fully  mixed 
throughout  the  soil  compartment  that  is  being  loaded  with  solid  phase  MC 
residue.  Thus,  a  continuously  stirred  tank  reactor  (CSTR)  model  is  assumed 
for  the  soil  model.  The  constituent  is  assumed  to  exist  in  solid  (free  product) 
and  non-solid  (dissolved  from  solid)  phases.  The  non-solid  phase  mass 
exists  in  equilibrium  distributed  as  dissolved  in  water  within  the  water  filled 
soil  pore  spaces,  as  adsorbed  from  water  to  soil  particles,  and  as  a  vapor  in 
air  within  the  air  filled  pore  spaces.  A  time-varying  mass  balance  is 
performed  for  both  the  solid  and  non-solid  phases. 

Solid  Phase  Mass  Balance 

The  solid  phase  (free  product  without  precipitation  from  aqueous  phase) 
mass  balance  is  stated  as 

^  =  m-Fdh-F„+FIndp  (2) 


where  Ms  is  the  solid  phase  mass  (g),  t  is  time  (yr),  L(t)  is  time-varying 
solid  phase  MC  mass  loading  (g/yr),  Fdis  is  dissolution  flux  (g/yr),  Fes  is  the 
erosion  flux  of  solid  phase  particles  (g/yr),  and  Fpretip  is  the  precipitation 
flux  (g/yr),  due  to  dissolved  pore  water  concentration  exceeding  the  water 
solubility  limit.  Loading  depends  on  munitions  use  each  year  as  described 
in  Chapter  3.  The  dissolution  and  precipitation  fluxes  will  be  described 
later  in  this  chapter,  and  the  solid  phase  particle  erosion  flux  is  described 
in  Appendix  A. 

Non-Solid  Phase  Mass  Balance 

The  non-solid  phase  mass  balance  is  stated  as 
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where  Mns  is  non-solid  phase  mass  (g),  Fr  is  the  rain-induced  pore  water 
ejection  and  runoff  flux  (g/yr),  Fe  is  MC  flux  due  to  soil  erosion  (g/yr),  Fi  is 
leaching  flux  (g/yr),  Fdecay  is  degradation  flux  (g/yr),  and  Fvoi  is  volatilization 
flux  (g/yr).  The  total  (particulate,  dissolved,  and  vapor)  non-solid  phase 
constituent  concentration  within  the  soil  matrix  on  a  total  volume  basis,  Ctt 
(g/ms)  is 

c  M„,  M,  Md  +  Mp  +  Ma 

"  V  AZb  AZb 

where  V  is  the  surface  soil  compartment  volume  (m3)  of  the  AOI.  The 
terms  Md,  Mp,  and  Ma  are  the  mass  dissolved,  mass  adsorbed  to  soil 
particles,  and  mass  of  vapor  in  air  (g),  respectively.  Note  that  the  surficial 
soil  layer  thickness  Zb  (m)  and  site  area  A  (m2)  are  assumed  to  be  constant 
in  time.  Thus,  it  is  assumed  that  an  active  soil  layer  of  the  same  thickness 
is  reestablished,  although  there  is  soil  loss  due  to  erosion.  The  total  non- 
solid  concentration  can  be  expressed  as  the  sum  of  the  three  phase  media 
concentrations  corrected  for  media  volumes  or  mass, 

c„=e„c,+(t-e, jcg+Ptca  (5> 


where 

Ci  =  concentration  dissolved  in  pore  water  or  liquid  phase,  Md/Vw, 
g/ms 

Cg  =  vapor  or  gas  concentration  in  air,  Ma/Va,  g/m3 
Ca  =  concentration  adsorbed  to  soil  particles,  Mp/Msoii,  mg/kg 
0W  =  soil  volumetric  moisture  content  or  ratio  of  water  volume  to 
total  volume,  fraction;  6W  can’t  be  greater  than  soil  porosity 
(f>  =  soil  porosity  or  ratio  of  void  volume  to  total  volume,  fraction 
pb  =  soil  dry  bulk  density,  g/ml  or  kg/L 
Vw  =  volume  of  water  in  the  AOI  soil,  m3 
Va  =  volume  of  air  in  the  AOI  soil,  m3 
Msoii  =  mass  of  soil  in  the  AOI,  g 

The  vapor  and  particulate  concentrations  can  be  related  to  the  liquid 
concentration  through  phase  equilibrium  partitioning, 

Ca=KdCt 


(6) 
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Cg  =  KHCl 


(7) 


where  Kd  (L/kg)  is  the  distribution  coefficient  for  partitioning  a  constituent 
between  soil  particles  and  water,  and  Kh  is  the  dimensionless  Henry’s 
constant  for  partitioning  between  air  and  water.  Kh  is  computed  from 
Henry’s  law  constant  He  (atm9m3/g-mol), 


Kh 


He 

RuTa 


(8) 


where  Ru  is  the  universal  gas  constant  =  8.206  E-5  atm0m3/g-mol  °K,  and 
Ta  is  ambient  absolute  temperature  (°K),  or  soil-water  matrix  temperature 
in  this  case,  which  is  the  temperature  in  °C  plus  273.  Substitution  of 
Equations  6  and  7  into  Equation  5  results  in 


C, 


CH 


0WR 


(9) 


where  the  retardation  factor  R  is  defined  as 

D_,Jt-Ow)KH+PbKd 


(10) 


Flux  Terms 

The  discussions  of  flux  terms  are  divided  among  those  that  are  more  easily 
described  (leaching,  degradation,  volatilization,  and  soil  erosion)  and 
those  that  have  a  more  involved  description  and  derivation,  which  include 
rainfall-induced  pore  water  extraction  and  runoff  and  dissolution.  Thus, 
three  flux  term  sections  follow  below.  The  solid  phase  constituent  particle 
erosion  flux  is  discussed  in  Appendix  A,  because  the  computation  of  this 
flux  term  requires  special  considerations  and  discussion. 

Leaching,  Degradation,  Volatilization,  and  Soil  Erosion  Fluxes 

The  leaching  flux  is  computed  from 


Fl=qwACl 


(11) 
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where  qw  is  the  average  annual  Darcy  water  infiltration  rate  (m/yr). 
Although  the  soil  model  is  time-varying,  annual  average  values  are  used 
for  the  hydrologic  fluxes,  which  greatly  simplify  model  input  without 
compromising  the  long-term  fate  forecasts  due  to  the  linearity  in  flux 
process  descriptors.  Substitution  of  Equation  9  into  Equation  11  results  in 


F  _  QwACtt 

1  ewR 


(12) 


The  degradation  flux  F decay  is  computed  from 

Fdma  =  AZt(XiewCl+XapbCa)  (13) 

where  Ai  and  Aa  are  the  degradation  rates  (yr1)  for  the  liquid  (aqueous 
dissolved)  and  aqueous  adsorbed  concentrations,  respectively.  It  is  assumed 
that  the  vapor  phase  does  not  degrade,  but  it  can  volatilize—  or  diffuse  — 
from  the  soil  into  the  overlying  air.  Using  Equations  6  and  9,  Equation  13 
can  be  rewritten  as 


F decay  ~ 


A'i 


The  volatilization  flux  Fvoi  is  computed  from 

Fvoi  =  KvA{0-ew)Cg 


(14) 


(15) 


where  Kv  is  the  volatilization  rate  (m/yr),  or  vapor  escape  rate  from  soil  to 
the  overlying  air.  The  calculation  for  Kv  is  discussed  later  below  in  a 
separate  section  following  the  Flux  Terms  section.  Using  Equations  7  and  9, 
Equation  15  can  be  rewritten  as 

FM=K,A^-e^Ka^-  (16) 

UwK 


Equations  12, 14,  and  16  can  be  rewritten  as  follows 


F  —Sdn.  a  F  C 

rl  ~  Q  ^ rdp  '-'ft 


(17) 
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Fdecay=AZb(AIFdpCtt+\aFppCtt) 

Fvol  —  Kv  A  Fap  Ctt 


(18) 

(19) 


The  quantity 


F 

0. 


is  the  pore  water  infiltration  rate.  The  factors  FdP,  Fpp, 


and 


Fap  are  used  to  convert  from  total  concentration  on  a  total  volume  basis  to 
dissolved,  particulate,  and  air  concentrations,  respectively,  in  the  porous 
media  on  a  total  volume  basis.  These  factors  are  defined  as 


F  -1 
dp  Sw+^-ew)K„+ptKd  R 


(20) 


p  _ _ Pb  K,i _ 

pp  ew+(t-ew)KH+PbKd 


(21) 


F  (0-ow)kh 

ap  0w+{^-9w)KH+pbKd 


(22) 


It  is  pointed  out  that  the  formulation  for  FdP  here  is  different  from  the  one 
used  in  the  steady-state  soil  model  (Dortch  et  al.  2009). 

The  erosion  flux  is  computed  from 

Fe  =  E  ACtt  (23) 


where  E  is  the  average  annual  soil  erosion  rate  (m/yr).  It  can  be  shown 
that  Fe  includes  the  flux  of  chemical  adsorbed  to  eroded  soil  particles  and 
pore  water  chemical  that  is  within  the  eroded  soil  layer. 

Rainfall-Induced  Pore  Water  Ejection  and  Runoff  Flux 

Rain-induced  pore  water  ejection  and  runoff  is  used  to  estimate  runoff  flux 
Fr.  Chemicals  can  be  transferred  from  soil  pore  water  to  overland  runoff 
due  to  rainfall  impacting  the  soil  surface,  even  when  there  is  no  erosion. 
The  event -based  runoff  mass  removal  rate  of  pore  water  RedQ  (g/sec)  due 
to  rain-induced  ejection  can  be  computed  (Gao  et  al.  2004)  from 


(24) 
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where  er  (m/sec)  is  the  rate  of  soil  pore  water  ejection  during  a  rainfall 
event,  and  Ce  is  the  rainfall  event’s  time-averaged  soil  pore  water  chemical 

concentration  (g/m3)  in  the  soil  exchange  layer  adjacent  to  the  overland 
water.  The  soil  water  ejection  rate  er  is  defined  (Gao  et  al.  2004)  as 


alfi 

Pb 


(25) 


where  /  is  the  rainfall  intensity  (m/sec),  a  is  the  soil  detachability  (kg/L), 
and  <f)  is  the  saturated  water  content,  which  is  the  soil  porosity.  The 
instantaneous  soil  pore  water  chemical  concentration  in  the  exchange 
layer  during  a  rainfall  event  can  be  approximated  (Gao  et  al.  2004)  by 

Ce~C0exp(-/3t)  (26) 


where  C0  is  the  soil  pore  water  concentration  below  the  exchange  layer  and 
is  equal  to  FdPCtt/6w  or  Cz,  t  is  time  (sec),  and 


P 


eF 


dp 


dO 


(27) 


where  de  is  the  soil  exchange  layer  thickness  (m);  j 6  has  units  of  sec-1.  The 
event  time  average  of  Equation  26  is 


Caf  eptdt 


C  = 


(28) 


where  T is  the  time  averaging  interval,  which  is  the  event  duration  (sec). 
Substituting  the  definition  of  C0  and  Equations  25,  27,  and  28  into 
Equation  24  results  in 


(29) 

The  goal  is  to  be  able  to  apply  Equation  29  to  develop  an  average  annual 
rainfall  extraction  of  pore  water  with  runoff.  To  do  this,  an  average  or 
typical  rainfall  intensity  /  and  associated  duration  T  are  required  for  use 
in  Equations  25,  27,  and  29,  resulting  in  an  average  event  runoff  flux, 
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Rcr/y .  Multiplying  Re^  by  T  results  in  the  average  event  pore  water  mass 
removed  with  runoff.  The  number  of  such  events  occurring  within  a  year, 
N  (events/yr),  can  be  multiplied  times  Re^T  to  obtain  the  pore  water 
mass  removed  with  runoff  per  year  or  the  annual  runoff  export  Fr.  The 
above  statements  are  expressed  in  mathematical  form  as 

Fr=Ade(l-e-pf)CttN  (30) 

The  values  used  for  I ,  T  ,  and  N  should  satisfy  the  following  relation, 

Pr— IT  n  (3i) 

where  P,  is  the  average  annual  rainfall  (m/yr),  not  total  precipitation.  T  and 
/  drop  out  of  Equation  30  when  f3  is  multiplied  by  T  since  /?  has  rainfall 
intensity  in  it,  which  is  a  function  of  T  ,  Pr,  and  N  (see  Equation  31).  Thus, 
the  average  event  intensity  and  duration  do  not  affect  the  annual  export; 
rather  it  is  the  annual  rainfall  and  number  of  rainfall  events  that  are 
important  for  computing  annual  rainfall  extracted  pore  water  runoff  export. 
A  reasonable  approach  is  to  count  the  number  of  days  per  year  that  rainfall 
occurs  to  approximate  AT.  Equation  30  can  be  rewritten  as  shown  below 
given  the  above  discussion, 

Fr  =  Ade  (l  -  e  K )  CttN  (32) 


where, 


a<t>FdPpr 

PMN 


(33) 


During  rainfall  events,  the  upper  soil  is  water-saturated,  and  0W  equals  <ft, 
which  affects  the  calculation  of  FdP  and  k.  Thus,  Equation  33  becomes, 


aPr 

(  \ 

pan 

,  ^  +  PbKd  , 

(34) 


The  result  computed  by  Equations  32  and  34  is  fairly  sensitive  to  all  input 
parameters.  However,  the  only  two  parameters  that  are  not  easily 
determined  are  de  and  a,  so  typical  values  found  in  the  literature  like  those 
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reported  by  Gao  et  al.  (2004)  must  be  used.  It  is  expected  that  both  of 
these  parameters  are  affected  by  soil  texture,  land  use  and  cover,  and 
possibly  soil  chemistry.  It  is  emphasized  that  Pr  in  Equation  34  is  rainfall 
rather  than  total  precipitation;  thus,  snow  is  not  included,  and  N  is  the 
number  of  rainfall  days,  not  the  number  of  precipitation  days. 

Dissolution  Flux 

The  objective  was  to  develop  a  general  dissolution  formulation  that 
provides  estimates  of  solid  phase  mass  transferred  to  the  aqueous  phase 
on  an  annual  basis  using  average  annual  input  conditions  so  that  long¬ 
term  (decadal)  forecasts  can  be  easily  conducted.  The  formulation  needed 
to  take  into  account  dissolution  for  the  total  solid  phase  mass  deposited 
within  the  AOI,  not  merely  an  individual  particle,  and  should  be  applicable 
to  various  types  of  solid  phase  MC,  such  as  metals  and  HE. 


Dissolution  of  a  solid  particle  in  water  can  be  described  as  a  diffusion 
process  (Cussler  1997)  driven  by  the  concentration  gradient  around  a  solid 
particle,  which  is  expressed  as 


dmi 

~df 


D. 


A -(c.-cj 


(35) 


where 

m;  =  solid  phase  mass  for  particle  i,  g 

Dw  =  chemical’s  aqueous  phase  diffusion  coefficient,  cm2/sec 
Asi  =  surface  area  of  the  individual  solid  phase  particle,  cm2 
Cs  =  aqueous  solubility  of  the  chemical,  g/cnD 
Cia  =  chemical  concentration  in  water  adjacent  to  the  particle,  g/cm3 
h  =  the  boundary  layer  film  thickness  around  the  particle,  cm 
t  =  time,  sec 

During  quiescent  conditions  of  no  water  movement,  such  as  with  no 
rainfall,  Cia  approaches  Cs,  and  there  is  no  dissolution.  Thus,  it  can  be 
assumed  that  dissolution  occurs  only  when  water  is  moving  past  the 
particles  as  when  it  is  raining,  which  is  the  case  with  the  drop  impingent 
model  developed  by  Taylor  et  al.  (2009a).  If  water  is  continually  replaced 
around  the  particle  when  it  is  raining  such  that  the  adjacent  water 
concentration  is  dilute,  Cia  can  be  assumed  to  be  much  smaller  than  Cs. 
Equation  35  simplifies  to 
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dmi 

~dt 


<CsAsi 


(36) 


where  Dw/h  has  been  parameterized  with  Q,  the  dissolution  mass  transfer 
rate  or  solid  phase  surface  removal  rate  (cm/sec).  The  product  of  C, and  Cs 
is  the  mass  dissolution  flux  rate  /3m  (g/cm2-sec).  The  mass  dissolution  flux 
rate  j 6m  lumps  various  processes  to  account  for  weathering,  corrosion,  and 
dissolution  and  depends  on  constituent  solubility,  which  for  metals  can 
depend  on  a  number  of  site-specific  conditions,  such  as  soil  redox 
potential,  pH,  soil  chemistry,  and  other  factors. 

The  change  in  total  solid  phase  mass  m  over  time  (g/sec)  for  all  particles  is 

=  -&s  am  =  -pmam  =  -Fdis  (37) 


where  a  is  the  average  specific  surface  area  (cm2/g)  of  the  solid  phase 
mass  and  depends  on  the  distribution  of  the  size  and  shape  of  the  solid 
phase  particles  and  the  constituent  solid  phase  density.  The  product  a  m  is 
the  total  surface  area  of  solid  phase  particles.  The  total  solid  phase  mass  m 
is  the  same  as  the  previously  defined  variable  Ms.  The  value  of  a  generally 
increases  as  the  particles  dissolve  and  total  mass  decreases. 

Previously  reported  dissolution  models  for  high  explosives,  such  as  those 
reported  by  Phelan  et  al  (2004),  Taylor  et  al.  (2009a),  and  Lynch  et  al. 
(2002),  can  be  recast  in  a  form  similar  to  Equation  37.  The  Phelan  model 
is  similar  to  Equation  35  with  the  mass  transfer  rate  Dw/h  determined 
from  soil  column  leaching  studies. 

The  Lynch  model  is 


=  -amPmeeT™  (38) 

at 

where  Tw  is  ambient  water  temperature  in  degrees  C,  and  6  here  is  an 
empirical  coefficient  to  account  for  temperature  effects  on  dissolution  rate 
due  to  changes  in  solubility.  Equation  38  is  already  in  the  same  form  as 
Equation  37.  The  values  of  6  and  the  mass  dissolution  flux  rate  /3m  depend 
on  the  explosive  constituent  and  were  determined  experimentally  from 
batch-mixing  studies.  Values  of  6  were  found  to  vary  between  0.0574  and 
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0.0903  per  °C,  and  /3m  was  found  to  vary  between  1E-6  and  7E-5  mg/cm2- 
sec.  The  Lynch  model  tends  to  produce  rapid  dissolution  rates  since  the 
fitted  Pm  parameter  was  developed  from  continuous  batch-mixing  studies. 


The  Taylor-Lever  model  (Lever  et  al.  2005  and  Taylor  et  al.  2009a)  is  based 
on  individual  impingent  rain  drops  and  for  small  particles  is  stated  as 


dm  _  nVi 
dt  std 


(39) 


where  Vi  is  water  layer  volume  (cm3)  coating  a  spherical  explosive  particle 
and  is  a  function  of  the  water  thickness  h  coating  the  particle;  td  is  the  rain 
drop  arrival  interval  (sec).  This  model  is  valid  for  td  >  h2/Dw,  which  is 
satisfied  for  the  experimental  results.  Values  for  h  were  determined  by 
fitting  laboratory  experimental  results  from  water  dropping  on  individual 
particles,  and  the  value  was  found  to  be  on  the  order  of  0.1  mm  (Taylor  et 
al.  2009b). 

For  a  thin  water  layer  h  and  relatively  large  solid  particles  (e.g.,  on  the  order 
of  1  mm  or  greater),  Vi  can  be  approximated  as  the  product  ASih  for  an 
individual  particle,  and  the  ratio  h/td  is  similar  to  a  surface  renewal  (or 
removal)  rate,  C,.  Thus,  Equation  39  can  be  recast  in  the  form  of  Equation  36 
and  Equation  37. 

Equation  37  is  used  as  the  basis  for  the  dissolution  formulation  within  the 
TREECS  Tier  2  soil  model,  except  that  it  is  cast  in  a  form  that  can  be  used 
to  estimate  annual  rates.  Consider  precipitation  rate  as  the  driving  force 
for  the  mass  transfer  or  water  removal  rate,  C,.  Recent  results  of  Taylor  et 
al.  (2009b)  indicate  that  there  is  a  linear  relationship  between  dissolution 
rate  and  precipitation  rate,  and  they  demonstrate  that  it  is  possible  to  link 
annual  HE  dissolution  flux  to  annual  precipitation.  Thus,  Equation  37  can 
be  recast  by  replacing  C,  with  average  annual  total  precipitation  Pt  (cm/yr 
or  m/yr  depending  on  the  units  of  a  and  Cs)  resulting  in 

Fdls  =  PtaMsCs  (40) 


Taylor  et  al.  (2009b)  included  snowfall  in  equivalent  rainfall  amount  in 
their  calculations  for  comparison  with  experimental  results.  Thus,  total 
precipitation  is  used  in  Equation  40.  Appendix  B  presents  the  application 
and  comparison  of  Equation  40  to  experimental  data  of  Taylor  et  al. 
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(2009b)  and  comparison  with  results  of  the  linear,  annualized,  Taylor- 
Lever  model  (Taylor  et  al.  2009b).  The  results  from  Equation  40  compare 
favorably  with  both  experimental  results  and  the  linear  Lever-Taylor 
model  as  discussed  in  Appendix  B. 

All  other  variables  in  Equation  40  have  been  previously  defined.  However, 
relationships  are  required  for  estimating  a  and  Cs.  Considering  two  solid 
particle  shapes,  spherical  and  cylindrical,  a  can  be  derived  for  each  shape 
as  a  function  of  particle  average  diameter  di  (cm)  resulting  in 

2  4 

a(  cylinder)  — - 1 -  (41) 

PsJ  Psmd; 


a(  sphere)  = -  (42) 

PsA 

where  psm  is  the  solid  phase  constituent  mass  density  (g/cm3),  and  l  is  the 
length  (cm)  of  a  right  cylinder.  The  particle  diameter  di  and  the  specific 
surface  area  a  vary  with  time  as  mass  is  dissolved  from  the  particles.  For  a 
constant  particle  density,  it  can  be  shown  (Phelan  et  al.  2002)  that  average 
particle  diameter  di  can  be  related  to  total  particle  mass  Ms,  which  varies 
over  time,  with  the  result 


d,(t)  =  d,(t-  At) 


M,{t-At) 


(43) 


where  t  is  time,  and  At  is  the  time  step  in  the  solution  procedure.  Since  Ms  is 
being  calculated  each  time  step,  di  and  thus  a  can  be  computed  each  time 
step.  The  value  of  x  in  the  exponent  of  Equation  43  is  3  for  spherical 
particles  and  2  for  cylindrical  particles  of  constant  length  /.  If  solid  phase 
mass  is  decreasing  with  time,  the  average  particle  diameter  decreases  with 
time  and  vice  versa.  However,  the  average  particle  diameter  at  the  new  time 
step  is  not  allowed  to  exceed  the  initial  average  particle  diameter.  Thus,  if 
MC  mass  is  continually  being  loaded  into  the  AOI  at  a  rate  such  that  solid 
phase  mass  is  continually  increasing,  then  the  average  particle  diameter  will 
be  held  constant  at  the  initial  value.  Also,  the  average  particle  diameter 
should  not  be  allowed  to  decrease  below  a  minimum  particle  diameter  (such 
as  1.0  E-9  m)  to  avoid  division  by  zero  in  Equation  42.  In  reality,  older 
particles  are  shrinking  and  are  smaller  than  more  recently  loaded  and 
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newer  particles.  The  only  way  to  improve  upon  the  present  simplification 
would  be  to  track  the  fate  of  the  individual  MC  mass  classes  deposited  in  the 
AOI  each  year.  This  could  be  a  model  refinement  for  future  consideration. 

The  primary  reason  for  included  two  particle  shapes  is  to  investigate  the 
effect  of  shape  on  results.  Usually,  the  preferred  choice  would  be 
spherically  shaped  particles.  There  is  no  guidance  for  using  one  shape  over 
the  other. 

For  organic  MC,  the  water  solubility  Cs  is  a  chemical  property  that  varies 
with  ambient  water  temperature.  Phelan  et  al.  (2002)  measured  Cs  for  the 
explosives  RDX  and  TNT  for  varying  temperature  and  came  up  with  the 
relationships 


T 

w 

Cs  (TNT)  =  20.176  +  36.295e22  061 


Cs  (RDX)  = - - - - - 

s  0.0804  —  0.0194  In 


(44) 

(45) 


where  Tw  is  water  temperature  (deg  C),  or  soil-water  matrix  temperature  in 
this  case,  and  Cs  has  units  of  mg/L.  The  statistical  r2  values  for  Equation  44 
and  45  are  0.994  and  0.997,  respectively. 

Since  Equation  40  is  applied  for  annual  inputs  to  produce  annual  fluxes,  the 
value  of  Cs  should  be  an  annual  average  value.  A  reasonable  and  simplifying 
assumption  is  that  soil  water  temperature  from  precipitation  is  approxi¬ 
mately  equal  to  air  temperature.  In  fact,  mean  annual  soil  temperature  in 
temperate,  humid,  continental  climates  can  be  approximated  by  adding  1  °C 
to  the  mean  annual  air  temperature  reported  by  standard  meteorological 
stations  (NRCS  1993).  Thus,  average  annual  air  temperature  is  needed  to 
compute  solubility.  Daily  average  air  temperatures  can  be  averaged  over 
each  year  and  over  all  the  years  for  the  period  of  record  to  obtain  average 
annual  air  temperature. 

Water  solubility  values  exist  for  many  organic  chemical  constituents  for  a 
standard  temperature,  such  as  20  or  25  °C.  However,  water  solubility  of 
metals  is  much  more  variable  and  complex  and  depends  on  the  soil-water 
mineral  speciation  and  chemistry,  including  pH,  total  dissolved  solids, 
alkalinity,  and  potentially  other  factors.  Estimating  solubility  for  metals  is 
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beyond  the  scope  of  this  report,  but  it  does  need  to  be  addressed  within  the 
TREECS  project  and  is  addressed  in  a  follow-on  report  pertaining  to  the 
Tier  2  proof-of-concept  application  (Dortch  et  al.  2011). 

Volatilization  Rate 

Volatilization  from  soil  has  been  modeled  by  Jury  et  al.  (1983)  and  others. 
Their  one  dimensional  (vertical)  soil  models  display  contaminant 
concentrations  that  vary  over  soil  depth  as  a  result  of  diffusion  within  the 
soil,  volatilization  across  the  soil-air  interface,  and  other  processes,  such  as 
leaching  and  degradation.  A  simpler  approach  is  needed  here  since  the  soil 
is  treated  as  a  single  homogeneous  layer,  which  precludes  computing  time- 
varying  concentrations  that  vary  with  soil  depth.  However,  implementing  a 
simpler  approach  is  difficult  since  volatilization  mass  transfer  rate  varies 
with  time  as  contaminant  mass  is  depleted  near  the  soil  surface  due  to 
volatilization. 


The  volatilization  mass  transfer  rate  Kv  (m/yr)  is  required  to  compute  the 
volatilization  flux  as  shown  in  Equation  19.  The  volatilization  rate  is  a 
transport  coefficient  across  the  soil-air  boundary  that  can  be  approximated 
as  a  diffusion  coefficient  of  the  vapor  divided  by  the  diffusion  length  or 
thickness  for  mass  transfer.  It  is  assumed  that  the  diffusion  coefficient  of 
the  vapor  in  air  above  the  soil  surface  is  much  greater  than  the  diffusion 
coefficient  of  the  vapor  within  the  soil  matrix  due  to  the  porosity  and 
tortuous  pathways  between  solid  particles  in  the  soil.  Thus,  diffusion  within 
the  soil  layer  is  the  limiting  factor  for  volatilization.  The  effective  diffusion 
rate  coefficient  for  a  vapor  within  air  spaces  of  the  soil  matrix  can  be 
estimated  from  (Millington  and  Quirk  1961), 


it-ej0'3 


(46) 


where  DG  is  the  effective  diffusion  coefficient  (m2/day)  for  a  vapor  in  soil, 
and  Da"  (m2/day)  is  the  constituent  vapor  or  gas  diffusion  coefficient  in  air. 


Values  for  Df'  are  available  for  many  of  the  constituents  in  the  constituent 

databases  that  will  be  made  available  within  TREECS.  If  values  are  not 
available  within  a  database,  then  DGr Tor  the  constituent  of  interest  at 

20  deg  C  can  be  estimated  (Lyman  et  al.  1982)  from 
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DaGir  =  8.64  D2 


(47) 


where  D2  =  0.102  cm2/sec,  which  is  the  known  diffusion  coefficient  for 
carbon  disulfide  in  air;  M2  =  76  g/mole,  which  is  the  molecular  weight  of 
carbon  disulfide;  and  Mi  is  the  molecular  weight  MW  of  the  constituent  of 
interest. 

The  volatilization  rate  is  then, 


D, 


(48) 


Kv  =  365 


d, 


V 


where  dv  is  the  diffusion  layer  thickness  (m)  in  the  top  of  the  soil  layer.  This 
is  different  from  the  boundary  layer  thickness  immediately  above  the  soil- 
air  interface,  which  is  much  smaller  and  can  vary  with  wind  speed  and  other 
factors.  The  air  boundary  layer  thickness  has  a  value  on  the  order  of  1  mm 
(0.001  m)  as  calculated  from  experimental  data  (Spencer  et  al.  1988  and 
Hanna  et  al.  1982).  With  the  lower  effective  diffusion  coefficient  and  the 
larger  diffusion  thickness,  it  is  easy  to  see  why  volatilization  mass  transfer 
rate  in  soil  is  the  limiting  factor,  and  the  mass  transfer  rate  above  the  soil- 
air  interface  is  relatively  unimportant  when  dealing  with  a  single  layer  for 
the  soil  model.  Additionally,  the  present  model  does  not  consider  vertical 
migration  of  liquid  phase  constituent  due  to  evaporation,  which  can  also 
affect  the  volatilization  flux  (Jury  et  al.  1983). 

The  present  model,  including  Equations  19,  22,  46-48,  and  the  complete 
mass  balance  equations  were  applied  for  several  chemicals  which  had 
volatilization  fluxes  reported  in  the  literature  that  were  computed  with  the 
Jury  model  (Jury  et  al.  1983)  or  a  similar  model  of  that  type.  The  Jury 
model  is  considered  a  standard  for  determining  volatilization  from  soil.  By 
comparing  results  of  the  present  model  against  reported  flux  rates,  it  was 
possible  to  determine  an  appropriate  value  for  the  diffusion  thickness  in 
soil  dv.  These  comparisons  indicated  that  a  value  of  about  0.4  m  provided 
reasonable  estimates.  These  comparisons  are  discussed  below. 


During  the  development  of  the  CALTOX  model  (http://www.dtsc.ca.gov/Assessing- 
Risk/ctox  dwn.cfmh  the  Jury  model  was  applied  for  trichloroethylene  (TCE)  and 
hexachlorobenzene  (HCB)  to  aid  in  the  development  of  a  simpler 
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volatilization  routine  for  CALTOX.  The  application  and  the  results  are 
described  within  the  CALTOX  documentation  report  (Office  of  Scientific 
Affairs  1993!  (http://www.dtsc.ca.gov/AssessingRisk/ctox  dwn.cfml.  Given  the  initial  soil 
concentration  of  too  pmol/cm3  uniformly  distributed  throughout  the  soil  of 
depth  10  m,  the  Jury  model  indicated  that  the  volatilization  flux  from  the 
soil  varied  from  approximately  6o  and  0.3  pmol/cm2-day  at  year  zero  to 
approximately  0.5  and  0.02  pmol/cm2-day  at  year  10  for  TCE  and  HCB, 
respectively.  The  model  described  above  (with  dv  =  0.4  m  for  both 
chemicals)  gave  fluxes  of  86  and  0.026  pmol/cm2-day  at  year  zero  and 
3.9  and  0.026  pmol/cm2-day  at  year  10  for  TCE  and  HCB,  respectively. 
Although  the  present  model  does  not  give  exactly  the  same  result  as  the 
Jury  model  results,  it  does  give  approximately  the  same  order  of  magnitude, 
which  is  encouraging  given  the  simplicity  of  this  model.  The  soil  conditions 
used  for  the  Jury  model  were  not  reported;  therefore,  they  were  assumed  to 
be  the  following  for  the  present  model  application:  pb  =  1.48  kg/L;  (f>  =  0.44; 
Qw  =  0.175;  and  fraction  by  weight  of  soil  organic  carbon,  foc  =  0.01.  The 
chemical  properties  used  for  TCE  and  HCB,  respectively,  were:  He  =  9.69 
E-3  and  1.67  E-4  atm-m3/gmole;  molecular  weight  =  131.4  and  284.8; 
organic  carbon  -  water  partition  coefficient,  Koc  =  67.7  and  3380  L/kg;  and 
zero  degradation  rates.  The  comparisons  may  have  been  closer  if  the  soil 
conditions  used  for  the  Jury  model  had  been  known. 

Jury  et  al  (1990)  present  Jury  model  volatilization  results  for  various 
chemicals  in  soil  with  and  without  a  clean  soil  overburden.  Results  for 
benzene  in  sandy  soil  with  no  overburden  (zero  soil  cover)  after  1,000  days 
were  used  to  compare  with  the  present  model.  The  soil  conditions  for  the 
Jury  model  application  were:  pi,  =  1.59  kg/L;  </>  =  0.4;  0W  =  0.18;  and  foc  = 
0.0075.  The  chemical  properties  used  for  benzene  in  the  Jury  model  were: 
He  =  5.29E-3  atm-m3/gmole;  molecular  weight  =  78.1;  Koc  =  80  L/kg;  and 
degradation  half  life  of  365  days.  The  soil  layer  thickness  without  cover 
and  contaminated  with  benzene  was  1.0  m.  Again,  the  value  used  for  dv  in 
the  present  model  was  0.4  m.  The  present  model  projected  that  91  percent 
of  the  benzene  initial  mass  was  volatilized  after  1,000  days.  This  result 
compares  closely  with  the  approximately  89  percent  reported  by  Jury  et  al. 
(1990). 

Phelan  and  Webb  (1997)  applied  the  Jury  model  to  validate  a  model  they 
developed  of  similar  design  for  buried  mine  detection.  Thus,  the  Phelan 
model  is  a  Jury  type  of  model.  They  applied  their  model  for  TNT  to  evaluate 
how  long  it  would  take  for  TNT  volatilization  flux  to  be  emitted  at  the  soil 
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surface  and  the  level  of  concentrations  to  expect  over  time.  The  present 
model  was  compared  with  surface  flux  computed  from  the  Phelan  model 
after  the  surface  flux  had  reached  nearly  steady-state  for  a  constant  source 
loading  from  a  buried  mine.  Results  are  compared  for  scenario  l  of  the 
Phelan  model  study,  which  was  for  the  mid-continent  conditions.  The  soil 
conditions  for  the  Phelan  model  application  were:  pb  =  1.5  kg/L;  <f>  =  0.44; 
0W  =  0.25;  and  foe  was  not  specified  since  the  IQ  was  specified.  The  chemical 
properties  used  for  TNT  in  the  Phelan  model  were:  He  =  1.44E-8  atm-ms/g- 
mole;  molecular  weight  =  222;  Kd  =  6.8  L/kg;  and  degradation  half  life  of 
180  days.  For  the  Phelan  model,  the  mine  was  located  between  a  soil  depth 
of  5  and  15  cm.  The  soil  layer  thickness  contaminated  with  TNT  for  the 
present  model  was  set  to  0.15  m.  A  constant  TNT  source  in  the  soil  was 
input  in  the  Phelan  model  at  a  soil  depth  of  10  cm  at  a  rate  of 
8.6E-6  pg/cm2-day.  This  same  source  load  was  input  to  the  single  15-cm- 
thick  layer  of  the  present  model.  Again,  the  value  of  dv  used  in  the  present 
model  was  0.4  m.  The  fact  that  the  soil  layer  thickness  is  less  than  dv  did  not 
impede  application  of  the  present  model. 

The  nearly  steady-state  surface  volatilization  flux  after  one  year  computed 
by  the  Phelan  model  was  about  1.0E-11  pg/cm2-day.  The  nearly  steady- 
state  flux  computed  by  the  present  model  was  1.4E-10  pg/cm2-day.  The 
present  model  is  nearly  steady-state  from  the  beginning  of  the  simulation 
since  there  is  no  layer  of  clean  soil  above  the  contaminated  soil  like  there 
was  in  the  Phelan  modeling.  The  Phelan  model  results  had  not  quite 
reached  steady-state  after  one  year  and  were  still  increasing  slowly.  Also, 
the  Phelan  model  includes  the  effects  of  rainfall-evaporation  cycles  on 
migration  in  the  soil,  whereas  the  present  volatilization  model  ignores  the 
effects  of  evaporation  on  volatilization.  The  present  model  agrees  with  the 
Phelan  model  flux  results  within  an  order  of  magnitude,  which  is  not  bad 
considering  the  simplicity  of  the  present  model. 

It  is  noted  that  a  two  film  (or  dual  boundary)  model  was  considered,  where 
both  the  air  side  and  the  soil  side  boundary  layer  thicknesses  and  diffusion 
coefficients  were  considered  like  resistors  in  series.  However,  including 
the  air  side  boundary  had  only  a  very  minor  effect  on  the  results  given  the 
greater  resistance  associated  the  much  smaller  diffusion  coefficient  and 
much  larger  diffusion  layer  thickness  on  the  soil  side. 

The  relatively  simple  volatilization  model  proposed  above  provides 
reasonable  orders  of  magnitude  for  volatilization  flux  for  a  wide  range  of 
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chemicals  with  widely  ranging  Henry’s  constants  and  for  the  same  diffusion 
layer  thickness  of  0.4  m.  Overall,  the  fluxes  estimated  with  the  proposed 
model  are  a  little  greater  than  those  reported  in  the  literature.  However, 
volatilization  will  be  a  minor  process  relative  to  other  loss  pathways  for  MC 
since  the  Henry’s  constants  are  fairly  small  for  most  MC  of  primary  interest, 
such  HE.  The  user  has  the  option  in  TREECS  to  turn  off  the  volatilization 
pathway  by  re-setting  the  estimated  volatilization  rate  Kv  to  zero. 

The  soil  layer  thickness  Zb  will  be  one  of  the  inputs  for  the  time-varying 
soil  model  of  TREECS  Tier  2.  The  soil  contamination  thickness  can  vary 
widely;  therefore,  this  input  can  be  uncertain.  The  sensitivity  of  this  input 
variable  on  model  results  will  need  to  be  investigated  during  model 
testing,  but  initial  testing  indicates  that  it  does  slightly  affect  export  fluxes. 
A  value  of  about  0.5  m  for  Zb,  or  even  0.4  m  like  dv,  seems  reasonable.  The 
value  of  Zb  will  be  an  input  provided  by  the  user.  If  a  value  of  Zb  is  used 
that  is  less  than  the  suggested  value  of  dv  (i.e.,  0.4  m),  the  model  will  still 
run  without  problems.  The  values  used  for  dv  and  Zb  are  completely 
independent.  The  soil  model  user  interface  will  allow  the  user  to  change 
the  value  of  dv  from  the  default  value  of  0.4  m  if  desired  for  computing  Kv. 

Equation  Summary 

Without  Solubility  Limits 

The  above  defined  flux  terms  can  be  used  to  solve  the  non-solid  phase 
mass  balance  (Equation  3)  in  terms  of  the  unknown  Ctt, 

dC 

AZ„ =  F*  -  Ad,  (1  -  e-)NCa  -  EAC„ 

-jrAF-*C"  ~  AZ„  (A:F*  +  X,F„)C„  ~  K,AFC,  -  F„„,p 

(49) 


Collecting  and  cancelling  terms  in  the  above  equation  results  in 


dCtt 

Fdis 

K 

dt 

~  AZh 

zb 

—  +  F  -L. 

z,.  +  oz„  dp 


w  b 


(50) 
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Equation  50  and  Equation  2  constitute  a  system  of  two  coupled  ordinary 
differential  equations  that  can  be  solved  for  the  two  unknowns,  Ctt  and  Ms, 
which  vary  over  time.  The  coupling  terms  are  the  dissolution  and 
precipitation  fluxes.  The  precipitation  flux  term  in  Equations  2  and  50  is 
zero  when  solubility  is  not  limiting.  After  solving  the  concentration  Ctt,  all 
of  the  various  flux  terms  defined  above  can  be  calculated  for  each  time 
point  and  output  as  a  time  series  for  use  as  input  loads  to  other  models, 
such  as  leaching  flux  used  by  the  vadose  zone  model. 

Including  Solubility  Limits 

A  check  is  required  at  the  end  of  each  time  step  update  to  see  if  the 
computed  value  of  the  soil  pore  water  concentration  Ci  exceeds  the 
solubility  of  the  constituent  in  water  Cs  (g/m3  or  mg/L).  If  the  value  of  Ci 
computed  from  Equation  9  exceeds  Cs,  then  the  precipitation  flux  Fprecip  is 
computed  from 


F, . ,,  = 


AZb(C,-Cs) 


precip 


At 


(51) 


where  At  is  the  time  step  for  the  most  recent  time  update  of  the  solution.  If 
Ci  is  less  than  Cs,  then  Fprecip  is  zero  for  the  next  time  step  update. 

If  pore  water  MC  concentration  Ci  exceeds  solubility,  then  precipitation 
flux  computed  with  Equation  51  is  used  in  Equations  2  and  50  for  the 
solution  at  the  next  time  step. 

The  two  differential  equations  (Equations  2  and  50)  should  be  solved 
using  the  fourth-order-accurate  Runga-Kutta  time  integration  method 
with  an  option  to  use  either  a  constant  or  variable  time  step  that  is 
automatically  computed  to  maintain  stability.  Due  to  the  variable  time 
step  option,  it  is  recommended  that  the  solution  to  equation  51  be  lagged 
with  the  solution  of  Equations  2  and  50  as  described  above  rather  than 
trying  to  iterate  the  solution  with  those  two  equations,  which  could  lead  to 
convergence  problems. 
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5  Vadose  Zone  and  Aquifer  Models 

The  groundwater  below  the  AOI  soil  layer  will  be  simulated  with  two 
models,  an  unsaturated  (partially  saturated)  zone,  or  vadose  zone  model, 
and  a  saturated  zone,  or  aquifer  model.  As  constituent  mass  leaches  from 
the  upper  soil  horizon,  it  enters  deeper  soil  layers  that  are  partially 
saturated  with  water,  and  the  only  processes  are  vertical  percolation, 
adsorption,  and  degradation.  After  passing  through  the  vadose  zone,  the 
constituent  mass  reaches  the  top  of  the  groundwater  table  and  is  then 
transported  horizontally  with  the  groundwater  flow.  A  schematic  for  the 
groundwater  transport  modeling  is  shown  in  Figure  4. 

The  Multimedia  Environmental  Pollutant  Assessment  System  (MEPAS) 
(Buck  et  al.  1995)  consists  of  various  models  of  reduced  form  for  computing 
multimedia  fate  and  transport,  human  exposure  concentrations,  and  human 
receptor  doses  and  health  risks.  The  MEPAS  groundwater  models  consist  of 
time- varying  contaminant  fate/transport  models  of  the  vadose  zone  and 
aquifer.  The  MEPAS  version  5.0  groundwater  models  (Whelan  and 
Castleton  2006)  will  be  used  within  Tier  2  of  TREECS  to  compute  fluxes 
through  the  vadose  zone  and  aquifer  and  resulting  aquifer  concentrations. 
The  scientific  documentation  of  the  MEPAS  groundwater  models  (partially 
and  fully  saturated)  is  provided  by  Whelan  et  al.  (1996).  The  following  is  an 
overview  of  the  two  models. 


AOI 

Withdrawal  Well 

Soil  Layer 

Partially  Saturated 

Vadose  Zone 

MC  Plume 

Groundwater  Flow 
- ► 

Saturated  Zone, 

Aquifer 

Figure  4.  Site  schematic  of  transport  from  soil  to  groundwater. 
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Vadose  Zone  Model 

The  vadose  model  is  based  on  the  one-dimensional,  advection-dispersion 
equation, 

S£+(J5L)(3£)=(A.)(aV).Ac  (5: 

dt  Rj^  dzJ  dz2 

where, 

C  =  dissolved  concentration  (mg/L) 

Dz  =  dispersion  coefficient  in  the  z  direction  (cm2  s1) 

Rju  =  retardation  factor  for  unsaturated  zone  (dimensionless) 
t  =  time  (sec) 

w  =  pore  water  velocity  in  the  vertical  (Z)  direction  (cm  s_1) 
z  =  distance  along  the  vertical  coordinate  (cm) 

X  =  degradation  or  decay  rate  (sec1) 

The  retardation  for  the  unsaturated  zone  is  computed  from 


J?  =1  +  ^K, 

fu  d 


(53) 


where, 


pb  =  bulk  density  (g  cm  3) 

Kd  =  equilibrium  partition  (or  distribution)  coefficient  (mL  g1) 
0W  =  volumetric  moisture  content  (fraction) 


e 


w 


<t> 


g  id 

K 

s 


1/m 


(54) 


where, 

(j)  =  total  porosity  (fraction) 

Ks  =  saturated  hydraulic  conductivity  (cm  s1) 
K(dM ,)  =  unsaturated  hydraulic  conductivity  (cm  s_1) 
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m  =  coefficient  dependent  on  soil  type  where  m  =  2b  +  3 
b  =  soil  type  coefficient,  which  is  an  empirical  parameter  relating 
soil  matric  potential  and  moisture  content  and  is  dependent 
on  soil  texture  (Values  are  automatically  provided  from  a  table 
for  soil  textures  within  the  model  user  interface  and  are 
viewable  from  the  interface  help  file.) 

If  the  infiltration  rate  is  less  than  the  saturated  hydraulic  conductivity, 
K(0W)  is  set  to  the  infiltration  rate,  and  the  vertical  pore  water  flow  velocity 
w  is  equal  to  the  infiltration  rate  divided  by  6W.  If  the  infiltration  rate  is 
equal  to  or  greater  than  the  saturated  hydraulic  conductivity,  K(0W)  is  set 
to  the  saturated  hydraulic  conductivity  Ks,  and  w  is  equal  to  Ks  divided  by 
0W.  The  excess  infiltration  can  be  diverted  to  interflow  as  explained  in 
Chapter  8.  The  moisture  content  calculated  from  Equation  54  is  not 
allowed  to  be  less  than  the  specified  soil  field  capacity.  Thus,  the  moisture 
content  varies  between  the  field  capacity  and  saturation,  or  total  porosity. 

The  vertical  dispersion  is  computed  from 

Dz=azW  +  Dmd  (55) 


where, 


az  =  dispersivity  in  the  z  direction  (cm) 

Dmoi  =  molecular  diffusion  coefficient  for  solute  in  water  (cm2  s1) 

The  molecular  diffusion  is  assumed  to  be  negligible  compared  to  dispersion. 
The  dispersivity  is  generally  proportional  to  the  distance  from  the  source. 
Thus,  for  computing  flux  through  the  vadose  zone,  this  distance  is  the 
thickness  of  the  vadose  zone  layer.  The  constant  of  proportionality  (i.e.,  the 
factor  multiplied  times  the  vadose  zone  layer  thickness  to  get  az)  is  on  the 
order  of  0.01. 

The  vadose  zone  model  computes  the  concentration  entering  at  the  top  of 
the  vadose  layer  to  provide  the  boundary  condition  needed  for  the  analytical 
solution  of  Equation  52.  This  concentration  is  calculated  from  the 
constituent  infiltration  load  (g/yr)  entering  the  vadose  zone  from  above 
divided  by  the  percolation  flow  rate  (m3/yr).  The  percolation  flow  rate  is 
based  on  either  the  infiltration  rate  or  Ks,  whichever  is  smaller.  The 
infiltration  rate  for  the  vadose  zone  is  estimated  from  the  infiltration  flow 


ERDC/EL  TR-11-2 


36 


rate  (m3/yr)  passed  to  the  model  divided  by  the  flow  area,  which  is  obtained 
from  the  product  of  the  length  (L/)  and  width  (  Wf)  of  the  AOI,  which  are 
passed  to  the  model  from  the  soil  model.  Thus,  the  percolation  flow  rate  is 
either  the  infiltration  rate  or  the  saturated  hydraulic  conductivity  times  the 
flow  area. 

Aquifer  Model 

The  primary  function  of  the  aquifer  model  is  to  compute  the  MC 
concentrations  within  the  groundwater  plume.  A  Cartesian  coordinate 
system  is  used  for  the  plume  calculation  as  shown  in  Figure  5  along  with 
location  for  a  withdrawal  well  (receptor  location).  The  AOI  source  zone  is 
also  shown.  The  origin  of  the  coordinate  system  is  located  at  the  center  of 
the  AOI  vertical  projection.  The  well  is  located  x  distance  downstream 
from  the  origin  or  the  Y  axis,  y  distance  laterally  from  the  plume  centerline 
(X  axis),  and  z  distance  vertically  from  the  water  table  surface.  The  model 
computes  aquifer  concentrations  at  the  well  location  (x,  y,  z). 


The  MEPAS  aquifer  model  solves  the  time-varying,  one-dimensional 
advective,  three-dimensional  dispersive  mass  transport  equation  for  solute 
movement  through  a  porous  medium  with  a  unidirectional,  constant  or 
steady-state,  uniform  flow  velocity  and  with  first-order  degradation/decay 
and  equilibrium  sorption  partitioning.  Other  assumptions  include  the 
following: 

•  The  groundwater  environment  is  initially  free  of  contamination. 

•  All  transport  media  properties  are  homogeneous  and  isotropic. 

•  The  aquifer  is  of  finite,  constant  thickness  and  of  infinite  lateral 
extent. 
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•  Drawdown  effects  of  withdrawal  wells  and  other  transient  stresses  on 
the  aquifer  are  not  considered  by  the  semi-analytical  solutions. 

•  Flow  velocities  are  provided  by  the  user. 

•  Density  differences  between  a  contaminant  plume  and  the  natural 
groundwater  are  negligible. 


The  governing  equation  for  saturated  groundwater  transport  becomes 


dc 


u  „dC.  f  Dx  ,  Dy  ,  Dz  „ d2C 


^  +  (^X^)  =  (^)(^)  +  (^)(^m^)(„  2 

dt  Rfs  dx  Rfs  dx2  Rfs  dy  Rfs  dz2 


)  -  AC 


(56) 


where 


R„=l+^Kd 


(57) 


and 

u  =  pore  water  velocity  in  horizontal  (X)  direction  (cm/sec) 

Rfs  =  retardation  factor  in  the  saturated  zone  (dimensionless) 

D  =  dispersion  coefficients  in  the  x-,  y-,  and  z-directions 
(cm2/sec) 

(j)e  =  effective  porosity  (fraction) 

All  other  variables  have  been  previously  defined.  The  pore  water  velocity  is 
equal  to  the  Darcy  groundwater  flow  velocity  divided  by  the  effective 
porosity.  The  Darcy  flow  rate  (nD/yr)  at  a  specific  location  down-gradient  is 
computed  as  the  product  of  the  Darcy  velocity  times  the  vertical  and  lateral 
extent  of  the  plume  that  is  continuously  spreading  due  to  lateral  and  vertical 
dispersion  as  it  moves  downstream.  The  Darcy  flow  is  assumed  to  account 
for  the  recharge  flow;  hence,  the  recharge  flow  from  the  vadose  zone  is  not 
added  to  the  Darcy  flow. 

An  equation  analogous  to  Equation  55  is  used  to  estimate  Dx,  Dy,  and  Dz, 
except  that  u  rather  than  w  is  used,  and  there  are  dispersivities  for  the 
three  directions.  The  longitudinal  dispersivity,  ax,  is  the  product  of  a 
distance  and  a  factor.  The  longitudinal  distance  is  x,  the  distance  from  the 
source  to  the  well.  The  longitudinal  dispersivity  factor  is  on  the  order  of 
0.1.  The  lateral  or  transverse  dispersivity  is  estimated  as  0.33 ax.  The 
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vertical  dispersivity  is  estimated  as  0.0025%.  Effective  porosity  is  the  pore 
space  that  contributes  to  advective  flow  that  is  not  dead-end  space. 

A  combination  of  analytical  and  numerical  methods  is  used  to  solve 
Equation  56  for  a  variety  of  boundary  conditions,  including  a  source  influx 
specified  over  an  area,  such  as  an  AOI.  The  solution  scheme  also  accounts 
for  boundary  effects  of  aquifers  of  limited  vertical  thickness  such  as  perched 
aquifers  or  aquifer  water  tables  that  are  relatively  close  to  bedrock. 

It  is  noted  that  neither  the  vadose  nor  aquifer  model  allows  input  of  initial 
concentrations  due  the  legacy  design  of  these  models.  This  limitation  is 
not  a  problem  for  the  intended  use  of  TREECS,  which  is  for  a  source  of  MC 
contamination  loading  into  the  vadose  zone  from  the  soil. 
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6  Surface  Water  Models 

A  wide  variety  of  surface  water  types  may  be  encountered  near  military 
installations,  such  as  streams,  lakes,  ponds,  estuaries,  and  wetlands.  The 
goal  is  to  select  a  model  type  that  can  be  used  for  any  type  of  surface  water. 
Two  model  types  were  considered  for  the  Tier  l  modeling.  These  were  the 
time-varying,  one-dimensional,  longitudinal,  Contaminant  Model  for 
Streams  (CMS)  described  by  Fant  and  Dortch  (2007),  and  the  RECOVERY 
model  (Ruiz  and  Gerald  2001),  which  is  a  time-varying  contaminant 
model  for  standing  or  pooled  surface  water.  The  RECOVERY  model  was 
selected  for  Tier  1  since  it  provided  reasonable  results  for  both  streams  and 
lakes  while  restraining  the  complexity  as  desired  for  Tier  1.  However,  for 
Tier  2,  both  CMS  and  RECOVERY  will  be  provided  so  that  both  flowing 
(streams  and  rivers)  and  standing  water  (ponds,  lakes,  wetlands,  etc.)  can 
be  better  represented.  Each  model  is  briefly  described  below. 

RECOVERY 

A  complete  description  and  documentation  of  the  RECOVERY  model  are 
presented  by  Ruiz  and  Gerald  (2001).  A  brief  description  is  provided  here  so 
that  the  differences  between  RECOVERY  and  CMS  can  be  better  contrasted. 
RECOVERY  simulates  the  long-term,  time-varying  concentration  of 
contaminants  in  surface  water  and  bottom  sediments  for  both  dissolved  and 
particulate  contaminants.  A  schematic  of  how  the  water-sediment  system  is 
handled  in  RECOVERY  is  shown  in  Figure  6.  The  water  column  is  treated  as 
a  fully  mixed  volume,  or  a  continuously  stirred  tank  reactor  (CSTR).  The 
bottom  sediments  are  divided  into  two  types,  a  surficial  mixed  sediment 
layer  at  the  sediment-water  interface,  and  deep  sediments  below  the 
surficial  mixed  layer. 

This  treatment  results  in  three  mass  balance  equations  with  three 
unknowns,  which  apply  to  the  water  column,  mixed  sediment  layer,  and 
deep  sediment  layers.  Two  coupled  ordinary  differential  equations  are 
solved  for  the  total  (particulate  and  dissolved)  contaminant  concentration 
in  surface  water  (Cw)  and  in  the  mixed  sediment  layer  (Cm).  A  partial 
differential  equation  is  solved  for  total  concentrations  over  sediment  depth 
[Cs(z)]  for  the  deep  sediment  layers.  The  deep  sediment  extends  below  the 
depth  of  contamination  into  clean  sediment  so  that  a  zero  concentration 
gradient  boundary  condition  can  be  applied  at  the  bottom  of  the  sediment 
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column.  Each  mass  balance  equation  accounts  for  mass  fate  processes, 
such  as  sorption,  degradation,  etc.,  which  are  shown  in  Figure  7. 
Equilibrium,  reversible,  sorption  partitioning  is  assumed.  The  equations 
and  solution  schemes  are  documented  by  Ruiz  and  Gerald  (2001). 
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Figure  7.  Fate  processes  within  the  RECOVERY  model. 


Loading  boundary  conditions  include  inflowing  contaminant  mass  due  to 
export  from  the  soil  model,  which  includes  rainfall  extraction,  erosion,  and 
interflow  fluxes.  There  is  also  an  option  to  enter  a  user-specified  constant 
external  loading  of  contaminant  mass,  such  as  those  due  to  wastewater 
discharges,  air  deposition,  or  stream  background  loadings.  Although  the 
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model  state  dependent  variables  are  total  concentration,  the  model 
produces  output  for  total  and  dissolved  concentrations  in  the  water 
column  and  bed  via  equilibrium  partitioning. 


CMS 


This  section  provides  an  overview  of  the  CMS.  The  complete  model 
description  and  documentation  is  provided  by  Fant  and  Dortch  (2007). 
The  CMS  is  a  one-dimensional  (lD),  longitudinal,  time-varying  model  for 
contaminant  fate  and  transport  in  streams  and  rivers.  The  model  assumes 
steady,  uniform  flow;  thus,  there  is  no  hydraulic  or  hydrologic  routing 
involved.  Two  mass  balance  equations—  one  for  the  water  column  and  one 
for  the  bed—  are  solved  for  two  unknowns,  total  (particulate  and  dissolved) 
contaminant  concentration  in  the  water  column  and  in  the  sediment  bed. 
Only  one  layer  is  used  for  the  bed  in  this  model.  The  model  outputs  total 
and  dissolved  concentrations  in  the  water  column  and  bed  via  equilibrium 
partitioning. 


The  modeled  fate  processes  include  water  column  advection  and  diffusion, 
partition  of  contaminant  between  sediment  and  water  (for  water  column 
and  bed),  settling  of  particulate  mass  from  the  water  column  to  the  bed, 
resuspension  of  sediment  bed  mass  into  the  water  column,  burial  of 
sediment  bed  material,  diffusion  of  dissolved  mass  across  the  sediment- 
water  interface,  degradation/decay,  and  volatilization.  The  model  solves 
the  one-dimensional,  longitudinal  advection-diffusion-reaction  (mass 
balance)  equation  for  the  constituent  concentration  in  the  water  column, 
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where, 

cb  =  concentration  of  the  constituent  in  the  sediment  bed  (M/L3), 

total  mass  on  a  total  volume  basis 
cw  =  concentration  of  the  constituent  in  the  water  column  (M/L3), 

total  mass  on  a  total  volume  basis 
Dx  =  stream  longitudinal  diffusion  coefficient  (L2/T) 
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Fdh  =  factor  to  convert  to  dissolved  constituent  concentration  in  the 

sediment  bed  pore  water,  dimensionless 
Fdw  =  fraction  of  the  constituent  dissolved  in  the  water  column, 

dimensionless 

Fpw  =  fraction  of  the  constituent  in  particulate  form  in  the  water 

column,  dimensionless 
FI  =  hydraulic  depth  of  the  stream  (L) 

kdw  =  decay  rate  of  dissolved  constituent  in  the  water  column  (T1) 
kpw  =  decay  rate  of  particulate  constituent  in  the  water  column  (T1) 
kv  =  stream  volatilization  rate  of  the  constituent  (L/T) 
t  =  time  (T) 

U  =  stream  mean  velocity  (L/T) 

Vd  =  mass  transfer  rate  across  the  sediment-water  interface 

resulting  from  diffusion  of  the  dissolved  constituent  (L/T) 

Vr  =  sediment  resuspension  rate  (L/T) 

Vs  =  suspended  solids  settling  rate  (L/T) 
x  =  downstream  distance  along  stream  (L) 

The  stream  reach  is  sub-divided  into  computation  nodes  spaced  equal 
distance  apart  along  the  stream  longitudinal  axis,  x.  Equation  58  is  solved 
for  each  computation  node.  The  sediment  bed  is  represented  as  a  single 
homogeneous  active  layer  (i.e.,  a  CSTR)  underneath  each  water  column 
node.  The  mass  balance  equation  for  each  nodal  bed  layer  is 

=  ~{kdbFdt  +  kpbFpb  )c 'b  +  —  ( Fdwcw  —  Fdbcb ) 


vr  ^  vs  vb 


(59) 


where, 

Fpb  =  fraction  of  constituent  in  particulate  form  in  the  bed 

h  =  active  sediment  layer  thickness  (L) 
kdh  =  decay  rate  of  dissolved  constituent  in  bed  pore  water  (T1) 

kpb  =  decay  rate  of  particulate  constituent  in  the  bed  (T1) 

Vh  =  active  sediment  layer  burial  rate  (L/T) 
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All  other  variables  have  been  defined.  Longitudinal  exchange  between 
adjacent  bed  segments,  such  as  pore  water  flow  or  diffusion  along  the 
stream  flow  axis  within  the  bed,  is  assumed  to  be  negligible. 

The  factor  and  fractions  in  Equations  58  and  59  are  defined  as  follows, 


F  =  1 

*  1 +  KJS 

(60) 

F  -  K^S 

PW  1  +  KdwS 

(61) 

F  -  1 

1 t  +  KJl-^p, 

(62) 

,,  ■') . 

'*  , 

(63) 

where, 

Kdw  =  sediment-water  partitioning  coefficient  for  the  constituent  in 
the  water  column  (L3/M) 

Kds  =  sediment-water  partitioning  coefficient  for  the  constituent  in 
the  sediment  bed  (L3/M) 

S  =  TSS  concentration  (M/L3) 

<f)  =  porosity  of  the  sediment  bed,  fraction 
ps  =  dry  sediment  particle  density  (M/L3) 

Constituent  portioning  between  sediment  and  water  is  assumed  to  be  in 
equilibrium,  linear,  and  reversible.  The  model  can  only  be  applied  to  a 
single  reach  where  hydraulic  properties  (£7,  Dx,  and  H)  are  assumed  to  be 
constant  and  uniform. 

The  model  allows  either  direct  input  of  the  volatilization  rate  or  calculation 
of  it  based  upon  Henry’s  Law  and  Whitman’s  two-film  theory  (Chapra 
1997)  as  modified  for  the  influence  of  water  flow  and  wind  as  described  by 
Fant  and  Dortch  (2007).  Mass  transfer  rate  of  dissolved  constituents 
across  the  sediment-water  interface  is  calculated  based  upon  flow 
conditions,  molecular  diffusivity  and  other  factors  as  described  by  Fant 
and  Dortch  (2007). 
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A  single  active  sediment  layer  is  modeled.  The  active  bed  layer  is  assumed 
to  have  constant  properties,  thus,  the  thickness,  volume,  porosity,  bulk 
density,  and  solids  mass  are  constant  over  time  (i.e.,  steady-state)  and 
assumed  to  be  uniform  over  the  stream  reach.  Solids  and  contaminant 
mass  that  are  buried  to  deep  sediments  are  assumed  to  be  lost.  Performing 
a  steady-state  solids  balance  for  the  bed  layer  results  in  the  following 
relationship, 


n  =  M_v 

Pb 


(64) 


where  all  variables  have  been  previously  defined.  The  dry  bulk  density  of 
the  bed  pb  can  be  computed  from  (l  -  (f)  ps .  By  specifying  any  two  of  the 

sedimentation  rate  terms  in  Equation  64  along  with  the  bulk  density  and 
TSS  (S),  the  remaining  rate  term  can  be  solved.  The  three  sedimentation 
process  rates  of  Equation  64  are  applied  throughout  the  modeled  stream 
reach  along  with  a  constant  and  uniform  value  of  S. 

An  implicit,  finite  difference,  numerical  solution  scheme,  based  on  a 
variation  of  the  Crank-Nicolson  (CR)  method  (Chapra  1997),  is  used  to  solve 
the  partial  differential  equation  for  surface  water  contaminant  concen¬ 
tration  (Equation  59).  The  CR  scheme  is  unconditionally  stable  even  for 
large  time  steps  and  reduces  numerical  dissipation.  A  choice  between  two 
numerical  integration  schemes,  Euler  and  fourth  order  Runge-Kutta,  is 
provided  to  solve  the  ordinary  differential  equation  for  sediment  bed 
contaminant  concentration  (Equation  59).  An  adaptive  time-stepping 
solution  option  has  been  implemented  for  solving  the  ordinary  differential 
equation  for  the  bed  to  ensure  stable  results. 

The  CMS  has  been  verified  against  known  analytical  solutions  (Fant  and 
Dortch  2006)  and  the  RECOVERY  model  (Fant  and  Dortch  2007).  The 
main  reason  for  using  CMS  is  for  cases  where  there  may  be  sediment  and 
water  contaminant  concentrations  that  vary  along  the  stream  reach  and 
over  time.  Although  the  CMS  was  developed  for  long-term  fate  of 
contaminants  in  sediment  and  water,  it  can  also  be  used  to  assess  short¬ 
term  spill  events  or  transient,  storm-water  transport.  The  CMS  is  also  better 
suited  than  the  RECOVERY  model  for  streams  that  are  tidally  influenced 
since  it  is  lD  and  includes  diffusion/dispersion,  which  can  be  used  to 
represent  tidal  mixing  in  lD  models  through  tidal  dispersion  (Thomann 
1972,  Fischer  et  al.  1979). 
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7  Model  Implementations 

The  Tier  2  soil  model  will  generate  outputs  (export  mass  fluxes)  that  are 
used  by  the  vadose  zone  model  and  the  surface  water  model  (RECOVERY 
or  CMS).  The  vadose  zone  model  will  generate  output  that  is  used  by  the 
aquifer  model.  If  groundwater  discharge  is  included,  then  the  aquifer 
model  will  generate  output  that  is  used  by  the  surface  water  model. 

This  chapter  discusses  the  input  requirements  for  each  model  and  how 
those  are  addressed  for  Tier  2  analyses.  Model  outputs  are  also  discussed. 
No  additional  assumptions  or  simplifications  are  made  for  the  receiving 
water  models  (vadose,  aquifer,  and  surface  water)  in  Tier  2.  Thus,  the  user 
will  have  full  control  to  use  the  entire  capabilities  of  each  receiving  water 
model. 

Soil  Model 

The  soil  model  described  in  Chapter  4  has  the  input  requirements  listed  in 
Table  1.  A  user  interface  will  be  developed  for  the  Tier  2  soil  model  for 
setting  up,  reviewing,  and  editing  model  inputs. 

There  will  also  be  several  options  in  the  soil  model  user  interface  that  the 
user  will  be  able  to  choose  when  applying  the  soil  model.  One  option  will 
be  to  include  or  not  include  solid  phase  particle  erosion  (i.e.,  Equation  Ay 
for  the  flux  term  in  Equation  2).  Within  the  user  interface,  the  user  can  use 
solubility  values  that  are  passed  from  the  constituent  database  or  input 
user  specified  values.  For  RDX  and  TNT  solubility,  there  will  also  be  the 
option  of  computing  it  from  Equations  44  and  45.  The  user  will  have  the 
option  to  calculate  or  input  a  user-specified  volatilization  rate  Kv  within 
the  user  interface.  If  the  calculate  option  is  chosen,  then  the  user  will  also 
have  the  option  of  changing  the  default  value  for  dv  (used  in  calculating  the 
volatilization  rate)  and  using  the  estimated  diffusion  rate  in  air,  a  database 
value,  or  a  specified  a  value. 

The  variables  listed  in  Table  1  are  discussed  in  Chapter  4.  These  variables 
fall  into  five  categories  including:  site  characteristics,  soil  properties, 
hydrology,  fate/transport  parameters,  and  chemical-specific  properties.  The 
sources  of  data  and  methods  for  estimating  these  inputs  are  discussed 
below. 
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Table  1.  Tier  2  soil  model  inputs. 


Parameter 

Units 

Description 

Site  Characteristics 

Lr 

m 

AOI  dimension  that  is  parallel  to  the  groundwater  flow  (This 
variable  is  not  actually  used  by  the  soil  model,  but  it  will  be 
read  by  the  soil  model  user  interface  since  it  is  required  by 
and  will  be  passed  to  the  vadose  and  aquifer  models.) 

Wf 

m 

AOI  dimension  that  is  perpendicular  to  the  groundwater 
flow  (This  variable  is  not  actually  used  by  the  soil  model, 
but  it  will  be  read  by  the  soil  model  user  interface  since  it 
is  required  by  and  will  be  passed  to  the  vadose  and 
aquifer  models.) 

A 

m2 

AOI  surface  area 

Zb 

m 

Active  soil  layer  thickness 

Tsoil 

°c 

Average  annual  temperature  of  the  soil-water  matrix, 
which  is  used  for  all  temperature  dependent  calculations 

L 

g/yr 

MC  mass  residue  loading  versus  time 

Cs(O) 

mg/kg 

Initial  solid  phase  MC  concentration  in  soil  on  a  soil  mass 
basis  at  time  0,  which  must  be  input  for  each  MC  that  is 
modeled  (This  variable  will  be  converted  by  the  soil  model 
user  interface  to  solid  phase  MC  mass,  g,  in  soil  for  use  in 
the  soil  model.) 

Cns(O) 

mg/kg 

Initial  total  non-solid  phase  MC  concentration  in  soil  on  a 
soil  mass  basis  at  time  0,  which  must  be  input  for  each 

MC  that  is  modeled  (This  variable  will  be  converted  by  the 
soil  model  user  interface  to  non-solid  phase  MC 
concentration  in  soil  on  a  total  volume  basis,  g/m3,  for  use 
in  the  soil  model.) 

Soil  Properties 

0w 

fraction 

Volumetric  soil  moisture  content 

pb 

g/cm3 

Soil  dry  bulk  density 

<]> 

fraction 

Soil  porosity 

Hydrology 

P 

m/yr 

Average  annual  precipitation 

Pr 

m/yr 

Average  annual  rainfall 

E 

m/yr 

Average  annual  soil  erosion  rate 

C|w 

m/yr 

Average  annual  water  infiltration  rate  through  soil 
compartment;  rate  of  vertical  water  movement  through 
soil  after  ET  and  runoff  losses  (equal  to  groundwater 
recharge  for  no  interflow) 

Qr 

m/yr 

Average  annual  surface  water  runoff  rate 

Fir 

fraction 

Fraction  of  annual  water  infiltration  flow  rate  and  mass 
flux  that  goes  to  soil  interflow 

N 

yr-i 

Average  number  of  rainfall  events  per  year 
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Parameter 

Units 

Description 

Fate/Transport  Parameters 

Kd 

L/kg 

Soil-water  constituent  partition  coefficient,  which  must  be 
input  for  each  MC  that  is  modeled  (Options  for  estimating 
this  parameter  will  be  provided  in  the  user  interface.) 

de 

m 

Soil  exchange  layer  thickness  for  rainfall  ejection  of  pore 
water 

a 

Kg/L 

Soil  detachability  for  rainfall  ejection  of  pore  water 

h 

yr-l 

Decay/degradation  rate  of  liquid  (water)  phase 
constituent,  which  must  be  input  for  each  MC  that  is 
modeled 

A,a 

yr-1 

Decay/degradation  rate  of  adsorbed  (particulate)  phase 
constituent,  which  must  be  input  for  each  MC  that  is 
modeled 

di 

pm 

Initial  mean  diameter  of  solid  phase  constituent  residue 
particles,  which  must  be  input  for  each  MC  that  is 
modeled 

Kv 

m/yr 

Volatilization  rate,  which  must  be  input  for  each  MC  that  is 
modeled  (The  soil  model  user  interface  will  have  the 
option  to  compute  or  specify  this  parameter,  and  the  value 
will  be  passed  to  the  soil  model.  If  this  parameter  is 
computed,  then  the  diffusion  layer  thickness  dv  (m)  will 
also  be  required  as  input  in  the  user  interface.) 

dv 

m 

Diffusion  layer  thickness  in  the  upper  soil  horizon,  which  is 
required  for  computing  volatilization  rate  in  the  soil  model 
user  interface  (The  default  value  is  0.4  m). 

X 

dimensionless 

Particle  shape  index,  which  is  used  to  specify  for  the  soil 
model  whether  cylindrical  (index  =  1)  or  spherical  (index  = 

2)  solid  phase  particles  are  being  considered  for  each 
modeled  MC  (A  value  must  be  input  for  each  MC  that  is 
modeled.) 

Chemical-specific  Properties 

Cs 

mg/L 

Aqueous  solubility  limit,  which  must  be  input  for  each  MC 
that  is  modeled 

He 

atm»m3/g-mol 

Henry's  law  constant,  which  must  be  input  for  each  MC 
that  is  modeled  (This  property  will  be  converted  by  the 
user  interface  to  dimensionless  units  for  use  by  the  soil 
model.) 

MW 

g/mol 

Molecular  weight  (molar  mass  or  average  molecular 
mass),  which  must  be  input  for  each  MC  that  is  modeled 
(This  variable  is  not  required  by  the  soil  model,  but  it  will 
be  required  by  the  soil  model  user  interface  for  computing 
volatilization  rate,  which  is  required  by  the  model.) 

Psm 

g/cm3 

Solid  phase  constituent  mass  density,  which  must  be 
input  for  each  MC  that  is  modeled 
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Site  Characteristics 

The  AOI  dimensions  and  surface  area  can  be  decided  within  the  GIS  module 
of  TREECS  or  other  similar  external  means.  The  area  of  the  AOI  (A)  used  in 
the  soil  model  is  not  necessarily  equal  to  the  product  of  the  length  (L/)  and 
width  (IV/)  of  the  AOI.  Figure  8  shows  a  hypothetical  AOI  with  site 
dimensions.  Both  length  dimensions  are  required  for  the  vadose  and  aquifer 
models  and  are  passed  to  them.  The  AOI  width  should  be  perpendicular  to 
the  groundwater  flow  direction  and  should  extend  across  the  full  AOI  area. 

It  can  generally  be  assumed  that  groundwater  flow  and  runoff  flow  are  in 
roughly  the  same  direction.  The  AOI  length  is  perpendicular  to  the  width. 
The  AOI  area  A  is  used  by  the  soil  model  and  must  be  computed  using  GIS 
or  other  methods  if  the  AOI  is  an  irregularly  shaped  polygon. 


The  active  soil  layer  thickness  Zb  is  the  upper  soil  horizon  where  most  of 
the  MC  mass  residue  is  believed  to  reside.  In  general,  this  thickness  can  be 
on  the  order  of  a  meter  or  less.  A  default  value  for  Zb  of  0.4  m  is  suggested, 
since  this  is  the  recommended  soil  diffusion  layer  thickness  for  computing 
volatilization  rates  (See  Chapter  4). 

The  average  or  mean  annual  temperature  of  the  soil -water  matrix  of  the 
active  soil  layer  is  required  to  compute  the  dimensionless  Henry’s  constant 
for  ambient  conditions  and  solubility  when  a  function  relating  solubility 
and  temperature  is  available,  such  as  for  TNT  and  RDX.  Mean  annual  soil 
and  soil-water  temperature  are  assumed  to  be  the  same  and  can  be 
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estimated  by  adding  l  °C  to  the  mean  annual  air  temperature,  which  can 
be  computed  from  daily  average  air  temperatures  within  the  HGCT. 

MC  mass  residue  loadings  are  required  and  are  passed  from  the  MC  residue 
loading  model  discussed  in  Chapter  3.  Additionally,  the  user  must  input  the 
initial  soil  concentrations  for  both  solid  (pure  product)  and  non-solid  phase 
(weathered)  MC.  These  concentrations  should  be  on  a  soil  mass  basis 
(mg/kg),  which  is  usually  how  soil  concentrations  are  reported.  The  solid 
phase  mass  balance  equation  requires  mass  rather  than  concentration. 
Therefore,  the  soil  model  user  interface  must  covert  the  solid  phase  (pure 
product)  mass  concentration  into  mass  (g)  units  by  multiplying  the  input 
initial  pure  product  concentration  by  the  soil  dry  bulk  density,  soil  layer 
thickness,  and  AOI  surface  area.  The  non-solid  phase  (weathered)  mass 
balance  equation  requires  concentration  on  a  soil  volume  basis  (g/m:>);  thus, 
the  soil  model  user  interface  must  multiply  the  input  initial  weathered 
concentration  by  the  soil  dry  bulk  density.  In  many  cases,  the  user  will  input 
zero  for  both  initial  concentrations.  When  there  is  reason  to  input  non-zero 
values,  such  as  when  the  initial  MC  total  (pure  product  plus  weathered 
phase)  concentration  in  soil  has  been  measured,  then  the  user  must  decide 
how  to  split  the  concentration  between  pure  product  and  weathered  phases. 
Since  soil  and  pure  product  MC  particles  usually  have  different  particle 
densities,  a  gravimetric  separation  technique  may  be  required  to  measure 
the  pure  produce  and  weathered  MC  concentrations.  In  the  absence  of 
special  analyzing  techniques,  the  split  between  the  two  types  of  MC  may 
have  to  be  based  on  a  guess  or  best  estimate. 

Soil  Properties 

Site-specific  soil  properties,  0W,  pb ,  and  (f>,  can  be  estimated  based  upon  soil 
texture,  such  as  silty  loam,  sandy  clay  loam,  etc.  (see  Dortch  et  al.  2009). 
Bulk  density  is  related  to  dry  soil  particle  specific  weight  or  density  ps  and 
soil  porosity  according  to  the  relationship  pb  =  (l  -  <j>)  ps .  The  value  of  ps  is 

typically  about  2.65  g/cm3.  Soil  texture  tables  include  field  capacity  rather 
than  moisture  content,  but  field  capacity  can  be  used  as  an  estimate  for  9W  in 
partially  saturated  surface  soils  since  the  water  content  tends  to  fluctuate 
above  and  below  the  field  capacity  over  the  period  of  a  year.  The  mean 
particle  diameter  of  the  soil  can  also  be  estimated  from  soil  texture  if  needed 
later.  The  Hydro-Geo-Characteristics  toolkit  (HGCT)  within  TREECS 
includes  tools  for  estimating  soil  properties. 
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Hydrology 

The  methods  being  built  into  the  HGCT  for  estimating  the  average  annual 
hydrologic  variables  (P,  P,,  E,  qw,  qr,  and  N)  are  presented  by  Dortch  et  al. 
(2009  and  2010).  The  average  annual  infiltration  ( qw )  and  runoff  (<+)  rates 
are  multiplied  by  the  AOI  area  in  the  soil  model  interface  to  obtain  flow 
rates  needed  by  the  downstream  models.  It  is  noted  that  the  runoff  is  not 
actually  used  in  the  soil  model  computations,  but  the  runoff  flow  rate 
(m3/yr)  is  calculated  and  written  to  the  WFF  surface  water  output  file  as  this 
information  can  be  consumed  by  downstream  surface  water  models,  such  as 
CMS.  Infiltration  rate  (m/yr)  is  needed  by  the  soil  model  and  infiltration 
flow  rate  (m°-/yr)  is  required  by  the  vadose  zone  model;  therfore,  the 
infiltration  flow  rate  is  written  to  the  WFF  vadose  output  file. 

Rainfall  is  precipitation  that  occurs  above  o  °C.  The  hydrologic  methods  in 
the  HGCT  require  a  long-term  (20  years  or  more),  historical  record  of 
daily  precipitation  and  monthly  average  air  temperatures.  Monthly 
average  air  temperatures  can  be  developed  from  daily  average  values. 

The  fraction  of  infiltration  flow  and  mass  flux  that  goes  into  soil  interflow, 
which  is  diverted  to  surface  water,  can  be  specified  by  the  user  based  upon 
their  best  judgment,  or  it  can  be  computed  by  the  soil  model  user  interface 
if  a  value  for  the  saturated  hydraulic  conductivity  for  the  vadose  zone 
below  the  soil  layer  is  provided.  The  method  for  computing  the  interflow 
fraction  is  explained  in  Chapter  8.  Interflow  water  flow  rate  is  added  to  the 
runoff  flow  rate  and  subtracted  from  infiltration  flow  rate  for  output. 
Likewise,  the  interflow  mass  flux  is  subtracted  from  the  leaching  flux  and 
added  to  the  pore  water  extraction/runoff  flux. 

Fate/Transport  Parameters 

There  will  be  several  options  for  selecting  Kd.  An  estimator  will  be  built 
into  the  Tier  1  and  2  soil  models’  user  interfaces  for  estimating  Kd  for 
organic  constituents.  The  estimator  (Streile  et  al.  1996)  is  based  on  soil 
texture  (percentages  of  sand,  silt,  and  clay)  along  with  percent  organic 
matter  and  the  organic  carbon  to  water  partitioning  coefficient  (Koc), 
which  can  be  estimated  from  the  octonol-water  partitioning  coefficient 
(Row)  if  not  available.  The  estimator  is  stated  as, 


Kd  =  0.0001  Koc  (57.735 OM  +  2.0 clay  +  0.4  siZt  +.005 sand )  (65) 
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where  OM  is  the  percent  by  weight  of  organic  matter  in  the  soil.  If  Koc  is 
not  known  but  Kow  is  known,  Koc  can  be  estimated  with  Koc  =0.617  foc , 

where /oc  is  the  fraction  by  weight  of  organic  carbon  in  the  soil.  Soil  organic 
carbon  is  about  58  percent  of  soil  organic  matter;  thus,  foc  =  0.0058  OM  . 


There  also  will  be  a  look  up  table  capability  for  a  limited  number  of 
inorganic  constituents  that  provides  Kd  as  related  to  soil  pH,  percent  clay, 
percent  OM,  and  percent  iron  and  aluminum  oxyhydroxides  (Strenge  et  al. 
1989).  If  there  are  soil  Kd  values  in  the  constituent  database,  those  can  be 
selected  as  well.  User-specified  values  for  Kd  can  also  be  input. 

Default  values  will  be  provided  for  four  of  the  other  fate/transport 
parameters  with  the  option  for  the  user  to  change  each.  The  default  values 
will  be  0.005  m,  0.4  kg/L,  0.0  yr1,  and  0.0  yr1  for  de,  a,  Si,  and  8a, 
respectively.  The  values  for  de  and  a  are  in  the  range  of  values  reported  by 
Gao  et  al.  (2004).  Zero  degradation  rates  are  for  conservative  conditions  of 
no  degradation.  Help  files  will  be  available  for  estimating  degradation 
rates  in  soil  if  non-zero  values  are  used.  These  help  files  will  be  based  on 
information  from  various  reports  (Jenkins  et  al.  2003,  Miyares  and 
Jenkins  2000,  Mulherin  et  al.  2005,  Gerald  et  al.  2007,  and  Dortch  et  al. 
2005).  Additionally,  the  constituent  databases,  especially  the  Army  Range 
Constituent  Database  (ARCDB),  may  have  degradation  half-  lives  for  MC 
in  soil.  If  there  is  a  soil  half-life,  the  soil  model  user  interface  will  pull  that 
value  into  the  application,  but  the  user  can  change  the  value  if  he  or  she 
desires.  The  soil  model  uses  degradation  rates  (yr1),  but  the  soil  model 
interface  takes  the  input  as  degradation  half-life  (yr)  and  converts  to  a  rate 
0  693 

using  A  =  — - ,  where  ti/2  is  half-life  (see  Equation  66).  It  should  be 

^1/2 

recognized  that  much  of  the  information  on  degradation  rates  is  based  on 
laboratory  studies.  Some  sites  can  experience  lower  degradation  rates  than 
reported  from  laboratory  studies  (Dortch  et  al.  2007). 

Solid  phase  particle  sizes  can  vary  widely  for  various  types  of  munitions 
and  for  a  particular  MC.  Help  files  will  be  provided  regarding  solid  phase 
MC  particle  sizes  for  setting  di.  These  help  files  will  be  based  on 
information  from  various  reports  (Larson  et  al.  2005,  Pennington  et  al. 
2005,  and  Taylor  et  al.  2004). 

There  will  be  a  feature  within  the  soil  model  user  interface  allowing  the 
user  to  either  specify  a  value  or  estimate  a  value  for  the  volatilization  rate 
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for  each  MC  using  Equations  46  and  48.  Of  course  the  volatilization  rate 
for  inorganic  MC  should  be  set  to  zero.  There  won’t  be  any  gas  phase 
existing  in  the  soil  for  inorganic  MC  since  Henry’s  constant  is  zero;  thus, 
there  should  not  be  any  volatilization  even  if  the  rate  is  not  set  to  zero.  If 
the  volatilization  rate  is  computed,  then  the  MC  diffusivity  in  air  must  be 
specified.  If  a  value  of  the  MC  diffusivity  in  air  exists  in  the  constituent 
database,  that  value  will  be  provided;  and  if  not,  a  value  will  be  provided 
using  Equation  47.  In  either  case,  the  user  can  over-ride  the  provided 
value  and  enter  his  or  her  value.  Additionally,  the  user  will  be  able  to 
either  specify  the  soil  diffusion  layer  thickness  for  the  volatilization  rate 
calculation  or  use  the  default  value  of  0.4  m. 

The  particle  shape  index  is  used  to  specify  for  the  soil  model  whether 
cylindrical  (index  =  1)  or  spherical  (index  =  2)  solid  phase  particles  are 
being  considered  for  each  modeled  MC.  For  cylindrical  and  spherical 
particles,  Equations  41  and  42,  respectively,  are  used  by  the  soil  model  to 
compute  particle  average  specific  surface  area.  The  value  of  x  in  the 
exponent  of  Equation  43  is  2  for  cylindrical  and  3  for  spherical  particles. 

In  most  cases,  spherical  particles  should  be  used.  The  primary  reason  for 
providing  cylindrical  particles  is  for  investigating  the  effect  of  particle 
shape  on  computed  results.  However,  cylindrical  particles  may  be  better 
suited  for  modeling  metal  projectiles  on  SAFRs. 

Chemical-Specific  Properties 

Constituent  chemical-specific  properties  will  be  provided  for  the  most  part 
by  databases  that  will  be  available  within  TREECS.  There  will  be  four 
options  for  selection  of  a  chemical  database  for  use  in  an  application.  One 
option  is  to  use  the  Army  range  constituent  database  (Zakikhani  et  al.  2002 
and  Dortch  et  al.  2005).  The  Army  range  constituent  database  (ARCDB)  can 
include  multiple  values  for  chemical-specific  properties,  so  a  viewer/editor 
is  required  and  will  be  provided  in  TREECS  to  show  the  user  each  value  and 
will  require  selection  of  a  value  to  use  for  each  property  and  for  each  MC. 
The  older  version  of  the  ARCDB  has  been  recently  updated  during  the 
TREECS  project  to  fit  the  new  specifications  and  schema  outlined  by  Dortch 
et  al.  (2005).  Additionally,  a  search  for  additional  property  values  is  being 
conducted  for  priority  MC  to  help  ensure  that  values  will  be  available. 

The  Framework  for  Risk  Analysis  in  Multimedia  Environmental  Systems 
(FRAMES)  (Whelan  et  al.  1997)  database  (http://mepas.pnl.gov/FRAMESVl/mmede.stm) 
and  the  Risk  Assessment  Information  System  (RAIS)  constituent  database 
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f http://rais.orni.gov/index.shtmO  will  also  be  available  in  TREECS  for  specifying 
constituent  chemical-specific  properties.  The  fourth  option  will  be  to  use  a 
user-defined  database.  User-defined  constituent  databases  are  based  on  the 
FRAMES  constituent  database  format.  A  constituent  database  editor  will  be 
provided  within  TREECS  tools  for  developing  and  editing  user-defined 
constituent  databases. 

Database  values  for  solubility  of  organic  MC  are  fairly  reliable.  Solubility  of 
metals  is  highly  dependent  on  the  site  soil  geochemistry.  The  equilibrium 
chemistry  model  Visual  Minteq  is  being  applied  for  estimating  metal 
solubility  given  soil-water  ion  (such  as  carbonate,  sulfate,  chloride,  calcium, 
magnesium,  and  potassium)  concentrations  and  other  soil  chemistry 
information.  The  results  of  the  investigation  of  using  Visual  Minteq  will  be 
presented  in  another  follow-on  report  to  this  report. 

Output 

Output  of  the  Tier  2  soil  model  will  consist  of  time  series  of  each  of  the 
following  computed  variables  for  each  MC: 

•  Solid  phase  MC  mass  (g)  in  AOI  soil,  Ms 

•  Non-solid  phase  mass  (g)  in  AOI  soil,  Mns 

•  Non-solid  phase  MC  total  (dissolved,  adsorbed,  and  vapor) 
concentration  in  AOI  soil  on  a  total  soil  volume  basis  (g/m3),  Ctt 

•  Non-solid  phase  MC  dissolved  concentration  in  AOI  soil  on  a  water 
volume  basis  (g/m3),  Ci 

•  Total  MC  concentration  (solid  and  non-solid  phases)  in  AOI  soil  on  a 
soil  mass  basis  (mg/kg),  Cts ■  These  concentrations  should  be  provided 
via  a  Soil  Concentration  File  (SCF)  according  to  the  FRAMES 
specifications.  Additionally,  the  concentrations  Ci  should  be  provided 
via  a  SCF  according  to  the  FRAMES  specifications.  FRAMES  is  used  as 
the  underlying  modeling  framework  that  provides  the  mechanisms  for 
interfacing  disparate  models.  The  user  needs  to  understand  nothing 
more  than  this  about  FRAMES. 

•  Mass  fluxes  (g/yr)  including: 
o  Dissolution,  Fdis 

o  Solid  phase  erosion,  Fes 
o  Non-solid  phase  erosion,  Fe 
o  Leaching,  Fi 

o  Rainfall  induced  pore  water  ejection  and  runoff,  Fr 
o  Degradation,  Fdecay 
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o  Volatilization,  Fvoi 

•  Vadose  WFF  MC  mass  flux,  g/yr,  according  to  the  FRAMES 
specifications.  The  flux  to  vadose  will  consist  of  Fi  less  the  portion  lost 
to  interflow  flux. 

•  Surface  water  WFF  MC  mass  flux  (g/yr)  as  adsorbed  (particulate)  and 
dissolved  according  to  the  FRAMES  specifications.  Particulate  flux  will 
consist  of  a  combination  of  Fes  plus  a  portion  of  Fe,  and  dissolved  flux 
will  consist  of  Fr  plus  a  portion  of  Fe  plus  the  interflow  portion  of  Fu 

Additionally,  the  soil  model  output  will  include: 

•  Vadose  WFF  water  flux  (percolation  flow),  m3/yr,  according  to  the 
FRAMES  specifications.  Percolation  flow  rate  (water  flux)  is  computed 
from  qw(\-Fif)A. 

•  Surface  water  WFF  water  flux  (runoff  and  interflow),  nH/yr,  according 
to  the  FRAMES  specifications.  Surface  water  flow  rate  (water  flux)  is 
the  sum  of  the  runoff  flow  rate  and  interflow  flow  rate.  The  runoff  flow 
rate  is  the  product  of  the  average  annual  runoff  depth  and  the  AOI  area 
A.  The  interflow  flow  rate  is  the  product  of  the  average  annual 
infiltration  rate,  the  fraction  diverted  to  interflow,  and  the  AOI  area,  or 
qwFifA  ■ 

Vadose  Zone  Model 

The  input  parameters  of  the  MEPAS  vadose  zone  model  are  listed  in 
Table  2.  The  inputs  are  divided  into  categories  according  to  the  input 
screens  of  the  user  interface  as  shown  in  Figure  9.  The  parameter  names  as 
identified  within  the  model  are  shown  in  parenthesis  in  the  table.  The 
units  shown  in  the  table  are  the  standard  units  for  model  input.  However, 
the  MEPAS  models  allow  multiple  input  units,  and  the  units  that  are  input 
are  automatically  converted  to  standard  units  required  for  the  model. 

The  AOI  length  dimensions  have  been  discussed  previously  for  the  soil 
model.  These  soil  model  inputs  are  passed  to  the  vadose  zone  model.  The 
groundwater  recharge  flow  rate  Qw  is  simply  the  net  infiltration  rate,  or 
percolation  rate,  which  is  qw  (1  -  Fif),  times  the  AOI  surface  area,  all  of 
which  were  input  to  the  soil  model.  The  variable  Fif  is  the  fraction  of 
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Table  2.  MEPAS  vadose  zone  model  input  parameters. 


Parameter 

Units 

Description 

Inputs  passed  from  soil  model 

Lf 

m 

AOI  dimension  that  is  parallel  to  the  groundwater 
flow 

Wf 

m 

AOI  dimension  that  is  perpendicular  to  the 
groundwater  flow 

WFF  aquifer 
water  flux,  Qw 

m3/yr 

Water  flow  rate  due  to  net  infiltration,  or  percolation, 
from  soil  (groundwater  recharge  into  vadose  zone) 

WFF  aquifer 
mass  flux,  (1- 
Fif)Fi 

g/yr 

MC  mass  flux  versus  time  due  to  leaching  from  soil 
to  vadose  zone,  which  is  the  total  leaching  flux  Fi 
times  (1  -  Fif),  where  Fif  is  the  fraction  of  infiltration 
flow  that  goes  to  interflow 

Soil  Composition 

WP-SAND 

Percent 

Percentage  of  sand 

WP-SILT 

Percent 

Percentage  of  silt 

WP-CLAY 

Percent 

Percentage  of  clay 

WP-OMC 

Percent 

Percentage  of  organic  matter 

WP-IRON 

Percent 

Percentage  of  iron  and  aluminum 

Characteristics 

pH  (WP-PH) 

pH  units 

pH  of  pore  water 

4>  (wp- 

TOTPOR) 

fraction 

Total  porosity 

0f  (WP- 

FIELDC) 

Percent 

Field  Capacity 

Ks  (WP- 
CONDUC) 

cm/day 

Saturated  hydraulic  conductivity 

Zv  (WP-THICK) 

cm 

Thickness  of  the  vadose  zone  layer 

az  (WP-LDISP) 

cm 

Longitudinal  (vertical  direction)  dispersivity 

pb 

g/cm3 

Soil  dry  bulk  density 

Constituent  Properties 

Kd  (WA- 
SUBKD) 

ml/g 

Adsorption  (partition)  coefficient 

Cs  (WP- 
(R)SOL) 

mg/L 

Water  solubility  of  constituent 

Ti/2  (WP- 
GHALF) 

days 

Half-life  of  constituent  in  vadose  zone 
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Figure  9.  MEPAS  vadose  model  user  interface. 

infiltration  flow  that  goes  to  interflow.  The  MC  mass  flux  due  to  leaching 
from  soil  that  enters  the  vadose  zone  is  Fi  (l  -  Fif),  where  Fi  is  the  total  soil 
leaching  flux.  The  soil  model  wrapper  will  compute  the  water  and  mass 
fluxes  entering  the  vadose  zone  and  pass  them  the  vadose  zone  model  in 
the  form  of  an  vadose  water  flux  file  (WFF). 

Inputs  Passed  from  Soil  Model 
Soil  Composition 

Soil  composition  includes  soil  texture  (percentage  of  sand,  silt,  and  clay), 
percentage  of  organic  matter  (OM),  and  combined  percentage  of  iron  and 
aluminum.  The  five  percentages  must  total  too  percent.  The  MEPAS 
models  include  a  utility  called  Soil  Class  for  selecting  the  soil  texture.  The 
user  interface  for  this  utility  is  shown  in  Figure  to.  Given  the  selection  of  a 
texture,  the  percentages  of  sand,  silt,  and  clay  are  automatically  filled  into 
the  data  boxes  in  the  user  interface  screen.  The  soil  texture  is  used  to 
automatically  set  most  of  the  parameters  on  the  Characteristics  screen. 

Percent  of  soil  organic  matter  can  be  estimated  as  OM  - 1 .72  OC ,  where 

OC  is  the  percent  of  soil  organic  carbon  by  weight  (g  organic  carbon/g  soil 
x  too).  Thus,  OM  can  also  be  estimated  from  OM  =  172  foc ,  where  /oc  is  the 

fraction  by  weight  of  organic  carbon  in  the  soil.  Organic  matter  is  used  for 
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Figure  10.  MEPAS  model  utility  for  selecting  soil  texture. 


estimating  Kd.  Iron  and  aluminum  content  refers  to  oxyhydroxides,  which 
are  important  relative  to  sorption  partitioning  for  some  inorganic 
constituents.  There  is  a  data  table  with  lookup  feature  in  TREECS  that 
relates  Kd  (for  a  select  list  of  inorganic  constituents)  to  soil  pH  and  the 
summed  percentage  by  weight  of  clay,  organic  matter,  and  iron  +  aluminum 
content  in  the  soil.  This  feature  was  briefly  described  for  the  soil  model 
inputs.  If  only  organic  MC  are  of  interest,  the  percent  of  iron  +  aluminum 
can  be  set  to  zero. 

The  soil  type  coefficient  is  also  shown  on  the  composition  input  screen,  but 
this  value  is  automatically  filled  in  when  the  soil  texture  is  completed, 
although  the  user  can  change  it.  The  user  must  also  input  the  pH  of  the 
soil  pore  water  for  use  in  computing  Kd  for  inorganic  constituents  that  are 
in  the  inorganic  database  as  explained  above.  Thus,  the  value  of  pH  is  only 
used  for  estimating  Kd  of  inorganic  MC  that  are  in  the  lookup  table 
described  above. 

Characteristics 

Once  the  soil  composition  screen  is  complete  through  the  use  of  the  soil 
class  input  screen,  all  but  two  of  the  parameters  on  the  characteristics 
screen  are  automatically  filled  in  based  upon  tables  that  relate  soil 
properties  to  soil  texture  (Dortch  et  al.  2009).  The  user  can  change  any  of 
the  soil  properties  on  the  characteristics  screen. 
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The  thickness  of  the  vadose  zone  layer  must  be  input  by  the  user.  This 
thickness  should  be  the  depth  from  the  top  of  the  vadose  zone  to 
groundwater,  since  TREECS  Tier  2  limits  the  vadose  zone  to  a  single  layer. 
It  should  be  noted  that  the  user  can  select  the  units  of  input,  which  is  true 
of  most  all  of  the  MEPAS  model  inputs. 

The  user  must  also  input  the  longitudinal  dispersivity,  which  is  actually  in 
the  vertical  direction  for  the  vadose  zone  model.  This  dispersivity  should 
be  approximately  0.01  times  the  vadose  zone  layer  thickness.  The  MEPAS 
models,  as  well  as  the  other  models  within  TREECS,  have  help  files 
accessible  from  the  user  interface  to  help  in  providing  inputs. 

Constituent  Properties 

Like  the  soil  model,  there  are  several  options  for  entering  the  soil 
partitioning  coefficient  Kd.  One  option  is  to  simply  enter  a  user-specified 
value.  Another  option  is  to  use  the  value  in  the  constituent  database  if  a 
value  is  present  in  the  database.  Another  option  is  to  estimate  Kd  based 
upon  soil  composition.  This  option  uses  Equation  65  for  organic 
constituents  or  the  look  up  table  for  inorganic  constituents  as  explained 
above  and  as  used  by  the  soil  model. 

Water  solubility  will  be  passed  to  the  model  interface  from  the  constituent 
database  that  is  used  within  TREECS  if  a  value  is  available  in  the  database. 
If  not,  the  user  must  enter  the  solubility. 


The  user  must  enter  the  half-life  in  vadose  zone.  Any  values  of  half-life  in 
soil  may  be  appropriate  for  the  vadose  zone  as  well.  Half-life  is  related  to 
decay  rate  A  (time1)  according  to  the  formulation 


L 1/2 


0.693 

A 


(66) 


where  the  time  units  (e.g.,  days,  years,  etc.)  can  be  selected  by  the  user. 

For  conservative  forecasts,  the  half-life  in  the  vadose  zone  could  be  set  at  a 
large  number,  such  as  1.0E20  years,  which  would  essentially  be  no 
degradation. 
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Output 

The  output  for  the  MEPAS  vadose  zone  model  include  time  series  of 
Aquifer  WFF  water  flux  (m3/yr)  due  to  groundwater  recharge  and  Aquifer 
WFF  MC  mass  flux  out  (g/yr)  associated  with  percolation.  This  mass  flux 
is  output  for  each  MC.  These  files  are  viewable  via  a  suite  of  viewers 
developed  for  FRAMES. 

Aquifer  Model 

The  input  parameters  of  the  MEPAS  vadose  zone  model  are  listed  in 
Table  3.  The  inputs  are  divided  into  categories  according  to  the  input 
screens  of  the  user  interface  as  shown  in  Figure  11.  The  parameter  names  as 
identified  within  the  model  are  shown  in  parenthesis  in  the  table.  The  units 
shown  in  the  table  are  the  standard  units  for  model  input.  However,  the 
MEPAS  models  allow  multiple  input  units,  and  the  units  that  are  input  are 
automatically  converted  to  standard  units  required  for  the  model. 

Inputs  Passed  from  the  Vadose  Zone  Model 

The  variables  passed  from  the  vadose  model  are  the  same  as  those  passed 
from  the  soil  model  to  the  vadose  zone  model,  except  that  the  mass  fluxes 
have  been  modified  by  fate/transport  through  the  vadose  zone.  Retardation 
and  degradation  can  greatly  alter  and  attenuate  the  time  series  of  mass  flux 
from  soil  as  it  is  passed  through  the  vadose  zone. 

Composition 

Soil  composition  is  input  in  the  same  manner  as  for  the  vadose  zone  model 
where  the  utility  is  available  in  the  user  interface  for  selecting  soil  texture. 
The  groundwater  soil  composition  can  be  set  differently  from  the  surface 
soil  and  vadose  zone. 

Characteristics 

The  percent  of  constituent  flux  entering  the  aquifer  is  usually  set  to  too 
percent  for  conservatism  in  predictions.  Values  less  than  too  percent  are 
intended  to  represent  interflow  that  goes  to  surface  water  rather  than 
groundwater.  However,  recall  that  interflow  will  be  handled  by  the  soil 
model  in  TREECS  before  infiltration  flux  enters  the  vadose  zone.  Thus, 
interflow  will  have  already  been  taken  into  account,  and  the  percentage 
entering  the  aquifer  from  the  vadose  zone  should  be  set  to  too. 
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Table  3.  MEPAS  aquifer  model  input  parameters. 


Parameter  Name 

Units 

Description 

Inputs  passed  from  vadose  zone  model 

Lf 

m 

AOI  dimension  that  is  parallel  to  the  groundwater  flow 

Wf 

m 

AOI  dimension  that  is  perpendicular  to  the  groundwater  flow 

WFF  aquifer  water 
flux 

m3/yr 

Water  flow  rate  due  to  percolation  through  the  vadose  zone; 
i.e.,  groundwater  recharge 

WFF  aquifer  mass 
flux 

g/yr 

MC  mass  flux  versus  time  due  to  percolation  from  the  vadose 
zone  to  the  aquifer 

Composition 

WZ-SAND 

Percent 

Percentage  of  sand 

WZ-SILT 

Percent 

Percentage  of  silt 

WZ-CLAY 

Percent 

Percentage  of  clay 

WZ-OMC 

Percent 

Percentage  of  organic  matter 

WZ-IRON 

Percent 

Percentage  of  iron  and  aluminum 

Characteristics 

WZ-FRAC 

Percent 

Percentage  of  constituent  flux  entering  the  aquifer 

WZ-PH 

Dimensionless 

pH  of  the  pore  water 

f  WZ-TOTPOR 

Percent 

Total  porosity 

(j)e,  WZ-EFFPOR 

Percent 

Effective  porosity 

Vd,  WZ-PVELOC 

cm/day 

Darcy  velocity 

Za,  wz-thick 

cm 

Thickness  of  aquifer 

pb,  WZ-BULKD 

g/cm3 

Soil  dry  bulk  density 

Concentration  Locations 

x,  WZ-DIST 

cm 

Longitudinal  distance  to  well 

y,  WZ-YDIST 

cm 

Perpendicular  distance  from  plume  center-line  to  well 

z,  WZ-AQDEPTH 

cm 

Vertical  distance  below  water  table  to  well  intake 

ax,  WZ-LDISP 

cm 

Longitudinal  dispersivity 

ay,  WZ-TDISP 

cm 

Transverse  dispersivity 

az,  WZ-VDISP 

cm 

Vertical  dispersivity 

Flux  Locations 

x,  WZ-DIST 

cm 

Longitudinal  distance  to  flux  location 

ax,  WZ-LDISP 

cm 

Longitudinal  dispersivity 

ay,  WZ-TDISP 

cm 

Transverse  dispersivity 

az,  WZ-VDISP 

cm 

Vertical  dispersivity 

Constituent  Properties 

Kd,  WA-SUBKD 

ml/g 

Sorption  partitioning  coefficient 

Cs,  WZ-RSOL 

mg/L 

Water  solubility 

T1/2,  WZ-GHALF 

days 

Haif-life  of  constituent  in  groundwater 
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Figure  11.  MEPAS  aquifer  model  user  interface. 

Once  the  soil  texture  is  selected  in  the  soil  class  tool  or  is  entered  on  the 
composition  screen,  the  total  porosity  and  dry  bulk  density  are 
automatically  filled  in  for  the  characteristics  screen.  However,  the  user  can 
enter  user-specified  values  if  desired.  If  pH  is  entered  by  using  the  soil  class 
tool,  it  too  shows  up  on  the  characteristics  screen.  Effective  porosity  must  be 
entered  by  the  user  and  is  equal  to  or  less  than  the  total  porosity.  A  help  file 
will  provide  suggested  values  of  effective  porosity  as  related  to  soil  texture. 
As  effective  porosity  decreases,  the  pore  water  velocity  increases.  However, 
the  retardation  factor  also  increases,  which  can  cause  the  constituent  front 
velocity  to  decrease. 

The  Darcy  velocity  is  the  most  sensitive  input  parameter  affecting  output 
of  the  aquifer  model  (Dortch  et  al.  2009  and  2010).  The  Darcy  velocity 
must  be  input  by  the  user,  but  the  HGCT  within  TREECS  can  be  used  to 
help  in  estimating  it. 

The  thickness  of  the  aquifer  moderately  affects  aquifer  model  output 
(Dortch  et  al.  2010).  This  parameter  must  be  input  by  the  user.  Site 
characterization  data  should  contain  information  regarding  the  depth  to 
groundwater  and  the  aquifer  thickness  or  depth  to  bedrock. 
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Concentration  Locations 

Tier  l  aquifer  model  only  allows  a  single  receptor  well  location  for  each  AOI 
application.  However,  Tier  2  will  allow  multiple  receptor  well  locations.  The 
longitudinal,  lateral,  and  vertical  distances  from  the  AOI  to  each  receptor 
well  location  must  be  input  by  the  user.  The  three  distances  are  illustrated 
in  Figure  5.  The  longitudinal  distance  x  is  the  distance  from  the  center  of  the 
AOI  to  the  well  in  the  down-gradient  direction  of  groundwater  flow.  The 
lateral  distance  y  is  the  perpendicular  distance  from  plume  center-line  to 
the  well,  and  the  vertical  distance  z  is  the  depth  from  the  water  table  surface 
to  the  well  intake.  To  ensure  conservative  groundwater  concentration 
predictions,  zero  values  can  be  used  for  y  and  z. 

Longitudinal,  transverse,  and  vertical  dispersivity  are  also  required  input, 
but  the  model’s  user  interface  can  estimate  these  values  for  the  user  if  the 
user  clicks  the  Estimate  button  on  the  user  input  screen.  Usage  location 
(well)  distances  and  dispersivities  must  be  input  for  each  usage  location, 
or  well.  As  before,  the  user  can  select  the  units  of  input  for  distance  and 
dispersivity. 

There  is  a  Flux  Location  tab  and  input  screen  in  the  user  interface  for 
groundwater  discharge  to  surface  water.  There  will  only  be  one  usage 
location  in  TREECS  for  fluxes.  The  flux  location  is  where  the  groundwater 
discharges  to  surface  water.  The  user  must  enter  the  longitudinal 
(horizontal)  distance  from  the  center  of  the  AOI  to  the  surface  water  body 
where  it  is  believed  that  groundwater  flow  leaves  the  aquifer  and  enters 
surface  water.  It  is  possible  that  the  target  water  body  (where  MC 
concentrations  will  be  evaluated)  may  be  further  downstream  of  the 
location  of  groundwater  discharge  to  surface  water.  This  is  allowed,  but 
there  will  be  no  MC  routing,  losses  or  attenuation  between  the  point  of 
groundwater  discharge  and  the  target  water  body.  In  other  words,  the  MC 
mass  flux  leaving  the  aquifer  enters  directly  into  the  target  water  body 
even  when  the  target  water  body  is  downstream  of  the  groundwater 
discharge  influent  location.  Like  the  Concentration  Locations  (i.e.,  well 
locations)  input  screen,  the  three  dispersivities  must  be  entered,  but  they 
can  be  estimated  by  the  model  user  interface.  The  user  can  change  the 
units  for  distance  and  dispersivity. 
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Constituent  Properties 

These  constituent  property  inputs  are  the  same  as  those  discussed  under  the 
vadose  zone  model.  It  should  be  recognized  that  Kd  and  half-life  values  can 
be  different  for  groundwater  compared  with  surface  soil  and  the  vadose 
zone.  The  OM  content  of  aquifer  soils  is  usually  much  lower  than  that  of 
surface  soils  and  possibly  the  vadose  zone,  thus  potentially  lowering  the  Kd. 
Some  constituents  degrade  more  slowly  in  groundwater  than  in  soil.  Others, 
like  RDX,  degrade  more  rapidly  under  anaerobic  conditions,  and 
groundwater  can  be  anaerobic. 

Output 

Output  of  the  MEPAS  aquifer  model  consists  of  the  following  time  series  at 
specified  locations  that  are  available  via  the  FRAMES  viewers: 

•  Aquifer  WFF  water  flux  (m3/yr)  in  the  groundwater  flow  at  the  flux 
location  for  groundwater  discharge  to  surface  water  or  at  the  first 
receptor  well  if  there  is  no  groundwater  discharge  to  surface  water 

•  Aquifer  WFF  MC  mass  flux  (g/yr)  associated  with  groundwater  plume 
for  each  MC  at  the  flux  location  for  groundwater  discharge  to  surface 
water,  or  at  the  first  receptor  well  if  there  is  no  groundwater  discharge 
to  surface  water 

•  Aquifer  WCF  (water  concentration  file)  groundwater  concentration 
(mg/L)  for  each  MC  at  each  receptor  well 

Surface  Water  Model  -  RECOVERY 

The  input  parameters  required  by  the  RECOVERY  surface  water  model  are 
listed  in  Table  4.  The  inputs  are  divided  into  categories  according  to  the 
input  screens  of  the  user  interface.  There  are  six  input  screens  with  the 
first  one  shown  in  Figure  12.  In  most  cases,  the  parameter  name/symbols 
shown  in  the  table  are  the  symbols  used  in  the  documentation  report  (Ruiz 
and  Gerald  2001);  otherwise,  brief  names  were  created  for  the  table.  The 
units  shown  in  the  table  are  the  standard  units  for  model  input.  The  user 
must  adhere  to  these  units  since  RECOVERY  does  not  allow  for  input  unit 
conversions  like  the  MEPAS  models  do. 
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Table  4.  RECOVERY  surface  water  model  input  parameters. 


Parameter  Name 

Units 

Description 

Inputs  Passed  from  soil  Model  or  Plus-SG  Operator 

WFF  (surface  water)  water  flux 

m3/yr 

Water  flow  rate  from  AOI  soil  runoff,  interflow,  and  groundwater 
discharge  to  surface  water,  which  are  combined  via  the  Plus- 
SG  Operator  (This  flow  is  provided  but  not  actually  used  by  the 
RECOVERY  model.) 

WFF  (surface  water)  mass  flux 

g/yr 

Combined  mass  fluxes  exported  from  AOI  soil  due  to  rainfall 
ejected  pore  water  runoff  (Fr),  soil  erosion  (Fe),  and  solid  phase 
particle  erosion  (Fes),  and  interflow.  This  mass  flux  can  also 
include  groundwater  mass  flux  due  to  groundwater  discharge 
to  surface  water.  Particulate  and  dissolved  fluxes  are  included 
in  the  Surface  Water  WFF,  where  the  dissolved  fluxes  include 
surface  runoff  (including  interflow),  dissolved  portion  of  soil 
erosion  mass  flux,  and  aquifer  dissolved  fluxes  which  are 
combined  by  the  Plus-SG  Operator.  The  particulate  flux  is  due 
to  the  particulate  portion  of  the  soil  erosion  mass  flux  plus  the 
solid  phase  mass  erosion.  The  RECOVERY  model  combines 
particulate  and  dissolved  fluxes  into  a  total  mass  loading  for 
model  calculations. 

Surface  Water  Morphometry  and  Flydrology 

Sw 

mg/L 

Total  suspended  solids  concentration  in  the  water  column 

foe  water 

Fraction 

Weight  fraction  carbon  in  solids  in  water  column 

Aw 

m2 

Long  term  average  water  surface  area 

Water  depth,  Hw 

m 

Longterm  average  of  surface  water  mean  depth 

Q 

m3/yr 

Average  annual  water  flow  through  rate 

Residence  time,  xw 

yr 

Surface  water  residence  time 

Mixed  Sediment  Layer 

L 

m 

Total  sediment  bed  depth  to  be  modeled 

Z 

m 

Depth  of  mixed  sediment  layer 

Am 

m2 

Mixed  sediment  layer  surface  area 

(j)m 

Fraction 

Mixed  sediment  layer  porosity 

pP  mixed  layer 

g/cm3 

Mixed  sediment  layer  particle  density  or  specific  gravity 

foe  mixed  layer 

Fraction 

Mixed  sediment  layer  weight  fraction  carbon  in  solids 

Deep  Sediment  Layers 

Fraction 

Deep  sediment  porosity 

pP  deep  sediment 

g/cm3 

Deep  sediment  particle  density  or  specific  gravity 

foe  deep  sediment 

Fraction 

Deep  sediment  weight  fraction  carbon  in  solids 

System  Properties 

Uw 

m/sec 

Mean  wind  speed 

Enhance  diffusion 

cm2/sec 

Enhanced  diffusion  between  mixed  sediment  layer  and  deep 
sediment 

Enhanced  mixing  depth 

cm 

Enhanced  mixing  depth  between  mixed  sediment  layer  and 
deep  sediment 
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Parameter  Name 

Units 

Description 

Vs 

m/yr 

Suspended  solids  settling  rate 

Vb 

m/yr 

Deep  sediment  burial  rate 

Vr 

m/yr 

Mixed  layer  sediment  resuspension  rate 

Constituent  Properties 

CwO 

Pg/L 

Initial  contaminant  concentration  of  constituent  in  water 

w 

kg/yr 

Additional  constant  external  loading  rate  of  constituent 

CmO 

mg/kg 

Initial  constituent  concentration  in  mixed  sediment 

CsO 

mg/kg 

Initial  constituent  concentration  in  deep  sediment 

Dm 

cm2/sec 

Molecular  diffusivity  in  water  of  constituent 

He 

atm-m3/g-mole 

Henry’s  Law  Constant  of  constituent 

MW 

g-mole 

Molecular  weight  of  constituent 

Kow 

(mg/m3  octanol)/ 
(mg/m3  water) 

Constituent  octanol-water  partition  coefficient 

Kdw 

L/kg 

Partition  coefficient  for  the  water  column  for  inorganic 
constituents 

Kdm 

L/kg 

Partition  coefficient  for  the  mixed  sediment  pore  water  for 
inorganic  constituents 

Kds 

L/kg 

Partition  coefficient  for  the  deep  sediment  pore  water  for 
inorganic  constituents 

kw  dissolved 

1/yr 

Decay  rate  for  dissolved  contaminant  in  water 

km  dissolved 

Vyr 

Decay  rate  for  dissolved  contaminant  in  mixed  layer 

ks  dissolved 

1/yr 

Decay  rate  for  dissolved  contaminant  in  deep  sediment 

kw  particulate 

1/yr 

Decay  rate  for  particulate  contaminant  in  water 

km  particulate 

1/yr 

Decay  rate  for  particulate  contaminant  in  mixed  layer 

ks  particulate 

1/yr 

Decay  rate  for  particulate  contaminant  in  deep  sediment 

Model  Control  Parameters 

T 

years 

Total  period  of  simulation 

N  print  1 

dimensionless 

Number  of  time  steps  between  print  intervals  for  output 

N  print  2 

dimensionless 

Number  of  time  steps  between  print  intervals  for  sediment 
layer  output 

N  layers  print 

dimensionless 

Number  of  sediment  layers  to  plot  in  output 
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Figure  12.  RECOVERY  surface  water  model  user  interface,  first  screen  of  six. 

Inputs  Passed  from  the  Soil  Model  or  Plus-SG  Operator 

The  only  output  passed  from  the  soil  that  will  be  used  by  RECOVERY  is  the 
Water  Flux  File  (WFF)  for  Surface  Water  time  series  of  mass  flux  exported 
from  soil  to  surface  water  and  from  groundwater  discharge  to  surface  water. 
The  soil  mass  fluxes  are  passed  from  the  soil  model  in  the  form  of  the 
surface  water  WFF,  which  has  adsorbed  (particulate)  and  dissolved  fluxes. 
The  particulate  flux  is  the  combination  of  soil  MC  mass  losses  due  to  soil 
particle  erosion  (the  particulate  portion  of  erosion)  and  solid  phase 
constituent  particle  erosion.  The  dissolved  flux  from  soil  is  due  to  rainfall 
induced  ejection  of  pore  water  and  its  runoff,  the  dissolved  portion  of 
erosion  flux,  and  soil  interflow.  The  Plus-SG  Operator  consumes  and 
combines  surface  water  WFF  and  aquifer  WFF  to  produce  surface  water 
WFF  (for  both  particulate  and  dissolved  portions)  which  is  consumed  by 
RECOVERY.  Thus,  the  discharge  of  groundwater  dissolved  constituent  to 
surface  water  is  added  to  the  dissolved  mass  flux  of  the  soil  model  WFF.  The 
distinction  of  WFF  dissolved  and  particulate  flux  is  required  for  compliance 
with  the  FRAMES  WFF  specifications.  The  RECOVERY  model  combines 
the  particulate  and  dissolved  loadings  into  a  single  total  loading  for 
compliance  with  RECOVERY  model  specifications.  The  Plus-SG  Operator  is 
discussed  further  in  Chapter  8. 
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The  RECOVERY  model  requires  the  average  annual  water  flow  rate  Q 
through  the  water  body.  The  soil  model  provides  an  average  annual  runoff 
flow  rate  (including  interflow),  but  this  flow  rate  may  be  quite  different  from 
the  flow  moving  through  the  water  body  due  to  runoff  from  other  drainage 
basins  in  addition  to  drainage  from  the  AOI.  Thus,  it  is  better  for  the  user  to 
have  to  provide  the  value  for  Q  based  on  site-specific  circumstances.  If  the 
AOI  provides  the  only  drainage  comprising  Q,  then  the  user  can  estimate  Q 
by  multiplying  the  AOI  average  annual  runoff  rate  by  the  AOI  area  and 
adding  the  interflow  rate  to  this  flow.  The  interflow  rate  is  the  excess 
infiltration  times  the  AOI  area  as  discussed  in  Chapter  8.  The  model  does 
not  require  the  MC  concentration  associated  with  water  entering  the  water 
body  via  the  water  body  through-flow.  If  there  is  a  background  loading  of 
MC  due  to  through-flow,  then  that  loading  can  be  included  in  the  Addition 
Constant  External  Loadings  field  of  the  constituent  properties  screen  of  the 
RECOVERY  user  interface.  The  MC  load  from  the  AOI  automatically  enters 
the  water  body  via  the  surface  water  WFF  file  that  RECOVERY  consumes. 

Surface  Water  Morphometry  and  Hydrology 

The  total  suspended  solids  concentration  of  the  water  column  (Sw)  is  a 
required  input  although  model  results  are  relatively  insensitive  to  this  input. 
The  value  for  Sw  should  be  the  long-term  average  concentration.  The 
fraction  by  weight  of  organic  carbon  of  the  suspended  solids  (foc  water)  is  a 
required  input.  The  water  column  foc  is  the  ratio  of  water  column  particulate 
organic  carbon  (POC)  to  total  suspended  solids  (TSS  or  Sw).  TSS  is  the  sum 
of  inorganic  suspended  solids  (ISS)  and  volatile  suspended  solids  (VSS), 
and  POC  is  approximately  0.4  VSS.  Water  column  foc  values  can  vary  widely 
depending  on  the  sediment  delivery  from  the  watershed  and  the 
productivity  of  the  water  body.  Based  upon  a  brief  review  of  observed  POC 
and  TSS  data  for  Chesapeake  Bay1  where  extensive  research  results  are 
readily  available,  values  of  water  column  foc  varied  from  below  0.01  to  more 
than  0.15,  and  typical  values  were  roughly  0.05.  Freshwater  systems  should 
have  similar  values. 

Three  of  the  following  four  variables  must  be  entered  by  the  user:  surface 
area  of  the  water,  Aw;  mean  water  depth,  Hw;  flow  rate  through  the  water 
body  (i.e.,  flushing  flow  rate),  Q;  and  water  residence  time,  xw.  These  four 
variables  are  related  as  follows, 


1  Personal  communication  with  Dr.  Carl  Cerco,  Research  Hydrologist,  Environmental  Laboratory,  U.S. 
Army  Engineer  Research  and  Development  Center,  Vicksburg,  MS. 
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Surface  water  flow-through  or  flushing  rates  Q  can  be  estimated  from  the 
product  of  the  catchment  basin  area  that  drains  into  the  surface  water 
body  and  the  estimated  average  annual  runoff  depth.  Procedures  for 
calculating  average  annual  runoff  depth  are  available  within  the  HGCT. 
Also,  the  estimated  mean  flow  velocity  in  the  water  body  can  be  multiplied 
by  the  mean  cross-sectional  flow  area  to  obtain  Q.  The  water  body  flow 
through  rate  should  represent  the  combined  flows  from  AOI  runoff, 
groundwater  discharge,  plus  all  other  sources  of  runoff  not  associated  with 
the  firing  range  or  AOI. 

Sediment  Mixed  Layer 

The  definitions  of  L  and  z  are  diagramed  on  the  second  input  screen  as 
shown  in  Figure  13,  where  z  is  the  depth  of  the  mixed  layer,  and  L  is  the 
total  sediment  depth  being  modeled  for  contamination.  There  are  model 
layers  below  L,  but  those  layers  are  assumed  to  be  clean  or  have  no 
contamination  initially.  A  typical  value  for  z  is  about  0.07  m  with  a  range 
of  roughly  0.05  to  0.2  m.  A  reasonable  value  for  L  is  1.0  m,  although  the 
mixed  sediment  layer  concentrations,  which  are  important  ecologically, 
are  not  very  sensitive  to  the  depth  of  L. 

The  mixed  sediment  layer  surface  area  (Am),  porosity  (())m),  specific  gravity 
or  sediment  particle  density  (pp),  and  foc  must  be  entered  by  the  user.  A 
good  assumption  is  to  set  Am  equal  to  Aw.  A  reasonable  value  for  mixed  layer 
porosity  is  0.7,  although  higher  values  (e.g.,  0.9)  are  possible  and  generally 
result  in  lower  mixed  layer  contaminant  concentrations.  Thus,  0.7  is  a 
conservative  value.  A  good  value  of  specific  gravity  for  most  sediment  is 
2.65.  The  mixed  layer  foc  is  the  ratio  of  mixed  layer  POC  concentration  and 
dry  sediment  bulk  density,  and  it  is  generally  less  than  foc  of  the  water 
column.  A  typical  range  of  mixed  layer  foc  is  from  less  than  0.01  to  up  to 
about  0.10  or  0.15.  Values  for  New  Bedford  Harbor  varied  from  a  maximum 
of  0.134  in  the  inner  harbor  to  a  minimum  of  0.001  in  the  middle  and  outer 
harbor  (Abdelrhman  et  al.  1998). 

Deep  Sediment  Layers 

Similar  to  the  mixed  layer,  the  user  must  enter  on  the  third  input  screen 
(which  is  not  shown  here)  the  porosity,  specific  gravity,  and  foc  for  the  deep 
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Figure  13.  Second  input  screen  of  the  RECOVERY  user  interface. 

sediments.  The  main  difference  between  screens  2  and  3  is  that  sediment 
vertical  profiles  can  be  input  for  the  parameters  of  the  deep  sediments. 
Consequently,  a  constant  value  for  all  deep  sediment  layers  can  be  entered 
for  each,  or  values  can  be  entered  for  each  layer  of  the  deep  sediments  for 
each  of  the  three  parameters.  A  typical  value  for  porosity  in  deeper 
sediments  is  about  0.5.  A  good  value  for  specific  gravity  is  2.65  as  before. 
The  deep  sediment  foc  is  generally  less  than  the  mixed  layer  value. 

System  Properties 

System  properties  are  entered  on  the  fourth  input  screen,  which  is  not 
shown  here.  Model  results  are  not  very  sensitive  to  wind  speed  since  most 
range  constituents  are  not  very  volatile.  Average  wind  speeds  are  on  the 
order  of  6  m/sec  for  the  U.S.  (Dortch  et  al.  2009). 

The  two  enhanced  mixing  parameters  are  in  the  model  to  emulate  the 
effects  of  bioturbation.  However,  in  most  applications,  it  is  reasonable  to 
set  both  parameters  to  zero. 
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Two  of  the  three  sedimentation  rate  parameters  (vs,  Vb,  and  vr)  must  be 
entered,  and  the  third  parameter  is  set  to  o  and  will  be  computed  by  the 
model  user  interface  using  a  steady-state  solids  balance  (Ruiz  and  Gerald 
2001).  The  settling  rate  (vs)  is  usually  the  easiest  one  to  estimate  and  is 
generally  on  the  order  of  30  to  100  m/yr  for  fine  grain  sediment.  If  the 
resuspension  rate  is  not  estimated,  it  could  be  assumed  to  be  very  small 
(e.g.,  1E-10  m/yr  or  practically  zero),  and  the  burial  rate  can  be  set  to  o  and 
will  be  calculated.  Resuspension  is  usually  fairly  small  for  deeper,  slowly 
moving  waters  that  are  wind-sheltered,  such  as  small  lakes  and  ponds. 

Constituent  Properties 

Constituent  properties  must  be  entered  on  the  fifth  input  screen  for  each 
MC  being  modeled.  The  button  labeled  “Next  Compound”  at  the  top  of  the 
screen  should  be  clicked  to  advance  to  each  MC  for  inputs  (see  Figure  14). 


Figure  14.  Properties  input  screen  of  the  RECOVERY  user  interface. 
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The  first  five  inputs  are  actually  constituent  initial  and  loading  conditions 
rather  than  properties.  For  most  TREECS  applications,  these  four  inputs 
(Cwo,  W,  Cmo,  and  cs0)  will  be  set  to  zero,  which  are  the  default  values  in  the 
user  interface.  These  zero  inputs  mean  that  there  is  no  contamination  in  the 
water  or  sediment  of  the  water  body  initially,  and  there  are  no  constant 
external  loadings  (such  as  atmospheric  deposition,  background  load  from 
through-flow,  or  point  source  discharges).  If  there  is  an  initial  contami¬ 
nation  of  MC,  values  for  sediment  and  water  concentrations  should  be 
entered.  Sediment  concentrations  are  much  more  important  than  water 
column  since  the  water  column  concentrations  are  flushed  out  quickly, 
relative  to  the  sediment.  Likewise,  if  there  are  constant  background  loadings 
of  MC  or  other  constant  external  loadings,  such  as  atmospheric  deposition, 
those  can  be  entered  as  well.  There  is  also  the  time-varying  loading  due  to 
the  AOI  export,  and  these  loadings  are  included  as  they  will  be  automati¬ 
cally  passed  from  the  soil  model  to  the  surface  water  model  via  the  WFF  as 
discussed  in  the  section  on  Inputs  Passed  from  the  Soil  Model.  There  can 
also  be  a  loading  due  to  groundwater  discharge  to  surface  water,  in  which 
case  the  Plus-SG  operator  combines  the  soil  model  WFF  with  the  aquifer 
WFF  as  explained  in  Chapter  8. 

The  constituent  properties  Dm,  He,  MW,  and  Kow  will  be  passed  from  the 
constituent  database  if  values  are  in  the  database.  Inorganic  constituents, 
such  as  metals,  should  not  have  database  values  for  He  and  KoW.  A  zero 
value  for  Henry’s  constant  will  be  used  in  those  cases,  and  volatilization 
flux  will  be  zero.  The  Kow  value  is  used  by  the  RECOVERY  model  to 
compute  Koc  from  Koc  =  0.617 Kow,  and  Kd  =/oc  K„c  for  organic  MC. 

Sediment  and  water  partitioning  coefficients  will  have  to  be  entered  for 
inorganic  constituents.  These  can  be  entered  by  clicking  the  “View 
Calculated  Data”  button  on  the  Properties  screen  of  the  user  interface  as 
shown  in  Figure  14  and  Figure  15.  Three  Kd  values  must  be  entered  for  these 
cases,  one  each  for  water  column  (Kdw),  mixed  sediment  layer  (Kdm),  and 
deep  sediment  layers  (Kds).  Selecting  appropriate  Kd  values  for  metals  can 
be  difficult  since  they  are  affected  by  the  local  geochemistry.  In  general,  Kdw 
and  Kdm  should  be  about  the  same  and  can  be  much  lower  than  Kds  due  to 
the  presence  of  sulfides  in  deeper,  anaerobic  sediments  that  can  bind  with 
metal  cations  forming  insoluble  precipitates.  Values  of  Kd  used  for  lead  and 
copper  in  the  Fort  A.P.  Hill  application  of  TREECS  Tier  1  were  4,000  and 
600  L/Kg,  respectively  (Dortch  et  al.  2010). 
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Figure  15.  Calculated  data  and  partitioning  coefficients  input  screen  of  the 
RECOVERY  model  interface. 

Six  different  decay  or  degradation  rates  (yr1)  can  be  entered  for  each 
constituent.  These  six  rates  are  for  dissolved  and  particulate  phases  for  the 
three  media  (water  column,  mixed  sediment  layer,  and  deep  sediments). 

For  conservatism,  the  default  values  are  zero.  Decay  rates  can  be  computed 
from  half-lives  (see  Equation  66).  Decay  rates  are  not  well  known  for  many 
organic  constituents,  so  if  the  values  are  set  above  zero,  sensitivity  should  be 
conducted  to  obtain  a  range  of  expected  results.  Decay  rates  for  most 
inorganic  constituents,  such  as  metals,  should  be  zero. 

Model  Control  Parameters 

Model  control  parameters  are  entered  on  the  sixth  and  final  input  screen, 
which  is  not  shown  here.  The  total  period  of  simulation  should  be  set  long 
enough  for  the  sediments  to  respond  to  changes  in  export  loadings  from 
soil.  For  example,  if  the  operational  period  of  the  firing  range  is  50  years, 
then  the  surface  water  model  simulation  period  should  be  set  to  at  least 
too  years.  Sediment  contaminant  concentrations  will  decrease  over  decadal 
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time  scales  following  cessation  of  range  use.  The  number  of  time  steps 
between  print  intervals  should  be  set  about  20  or  more.  The  value  only 
affects  the  quality  of  the  output  concentration  versus  time  plots.  The  same  is 
true  for  N print  2  except  that  this  parameter  is  used  for  deep  sediment 
animations  (concentration  profiles  in  sediment  over  time)  rather  than 
concentration  versus  time  plots  like  N print  1.  The  number  of  layers  to  plot 
is  also  used  for  the  deep  sediment  animations.  Each  layer  is  1.0  cm,  so  if  the 
deep  sediment  of  interest  is  0.5  m  deep,  N  layers  print  should  be  set  to  50. 

Output 

The  output  of  the  RECOVERY  model  includes  the  following  for  each  MC, 
which  are  viewable  via  a  special  viewer  developed  for  RECOVERY: 

•  Water  column  total  (particulate  plus  dissolved)  concentration,  pg/L, 
versus  time 

•  Surficial  (mixed  layer)  sediment  total  concentration  on  a  sediment 
mass  basis,  mg/kg,  versus  time 

•  Surficial  sediment  pore  water  dissolved  concentration  on  a  water 
volume  basis,  mg/L,  versus  time 

•  Volatilization  flux  from  water,  pg/yr,  versus  time 

•  Dissolved  mass  flux  between  sediment  and  water  column,  pg/m2-yr, 
versus  time 

•  Time  animation  plots  of  deep  sediment  concentration  profiles  over 
depth  for 

o  Total  concentration  on  a  sediment  mass  basis,  mg/kg 
o  Total  concentration  on  a  sediment  volume  basis,  pg/nD  and  mg/L 
o  Dissolved  concentration  per  volume  of  pore  water,  mg/L 

In  addition  to  the  RECOVERY  viewer,  there  are  three  general  types  of 
FRAMES  viewers  that  will  provide  output  results.  These  include: 

•  SCF  (Sediment  Concentration  File)  graphical  view  -  time  series  of 
surficial  (mixed  layer)  sediment  total  concentration  on  a  sediment 
mass  basis  (mg/kg)  and  surficial  sediment  dissolved  concentration  on 
a  water  volume  basis  (pore  water  concentration,  mg/L) 

•  WCF  graphical  view  -  time  series  of  water  column  total  concentration 
(mg/L)  and  dissolved  concentration  (mg/L) 

•  WFF  graphical  view  -  time  series  of  water  column  particulate 
(adsorbed)  mass  flux  (g/yr)  and  dissolved  mass  flux  (g/yr)  associated 
with  flow  of  water  with  constituents  through  the  water  body 
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Surface  Water  Model  -  CMS 

Only  one  surface  water  model  can  be  used  for  each  AOI  application;  thus, 
the  model  will  be  either  RECOVERY  or  CMS,  depending  on  the  type  of 
water  body.  The  input  parameters  required  by  the  CMS  surface  water 
model  are  listed  in  Table  5.  The  inputs  are  divided  into  categories 
according  to  the  input  screens  of  the  user  interface.  There  are  five  input 
screens  with  the  first  one  shown  in  Figure  16.  However,  when  operating 
within  TREECS,  the  last  input  screen  for  Loading  Data  will  not  need  user 
inputs  since  loadings  will  be  automatically  passed  via  the  WFF  from  the 
soil  model.  In  most  cases,  the  parameter  name/symbols  shown  in  the  table 
are  the  names  and  symbols  used  in  the  documentation  report  (Fant  and 
Dortch  2007);  otherwise  brief  names  and  symbols  were  created  for  the 
table.  The  units  shown  in  the  table  are  the  standard  units  for  model  input. 
However,  like  the  MEPAS  models,  the  user  can  use  different  units,  and  the 
model  user  interface  automatically  converts  them  to  model-standard  units. 

Inputs  Passed  from  the  Soil  Model  or  Plus-SG  Operator 

Similar  to  RECOVERY,  CMS  consumes  a  surface  water  WFF.  However, 
unlike  RECOVERY,  CMS  uses  both  the  water  fluxes  and  mass  fluxes  in  the 
WFF,  whereas  RECOVERY  only  uses  the  mass  fluxes,  or  loads.  The  water 
fluxes  are  added  in  the  CMS  to  the  constant  background  flow  rate.  The 
model  also  allows  background  loadings  in  addition  to  the  loadings  imported 
via  the  WFF.  Like  RECOVERY,  CMS  combines  WFF  particulate  and 
dissolved  mass  loadings  into  a  single  total  loading  for  model  calculations.  If 
there  is  groundwater  discharge  to  surface  water  when  using  CMS,  then  the 
Plus-SG  Operator  consumes  and  combines  surface  water  WFF  and  aquifer 
WFF  to  produce  surface  water  WFF  which  is  consumed  by  CMS,  just  like 
what  is  done  for  RECOVERY. 

Model  Parameters 

The  user  has  the  option  of  inputting  either  the  number  of  computational 
segments  or  the  length  of  each  segment.  Whichever  is  used,  the  other 
variable  is  computed  for  use  in  the  model  based  upon  the  distance  to  the 
usage  (or  receptor)  location  of  interest.  Typically  it  is  good  to  have  at  least 
20  segments  in  a  study  reach. 
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Table  5.  CMS  surface  water  model  input  parameters. 


Parameter  Name 

Units 

Description 

Inputs  passed  from  soil  model  or  Plus-SG  Operator 

WFF  (surface  water) 
water  flux 

m3/yr 

Water  flow  rate  from  AOI  soil  runoff,  interflow,  and  groundwater 
discharge  to  surface  water,  which  are  combined  via  the  Plus-SG 
Operator  (This  flow  is  used  by  the  CMS.) 

WFF  (surface  water) 
mass  flux 

g/yr 

Combined  mass  fluxes  exported  from  AOI  soil  due  to  rainfall 
ejected  pore  water  runoff  (Fr),  soil  erosion  (Fe),  and  solid  phase 
particle  erosion  (Fes),  and  interflow.  This  mass  flux  can  also 
include  groundwater  mass  flux  due  to  groundwater  discharge  to 
surface  water.  Particulate  and  dissolved  fluxes  are  included  in 
the  Surface  Water  WFF,  where  the  dissolved  fluxes  include 
surface  runoff  (including  interflow),  dissolved  portion  of  soil 
erosion  mass  flux,  and  aquifer  dissolved  fluxes  which  are 
combined  by  the  Plus-SG  Operator.  The  particulate  flux  is  due  to 
the  particulate  portion  of  the  soil  erosion  mass  flux  plus  the  solid 
phase  mass  erosion.  The  CMS  combines  particulate  and 
dissolved  fluxes  into  a  total  mass  loading  for  model  calculations. 

Model  Parameters 

Number  of  segments 

NA 

Number  of  computational  segments  in  the  modeled  stream 
reach 

Segment  length 

m 

Length  of  each  computational  segment 

At 

day 

Model  time  step 

TP 

yr 

Total  period  of  time  for  the  simulation 

Dx 

m2/day 

Longitudinal  dispersion  coefficient 

TSS 

mg/L 

Total  suspended  solids  concentration  in  the  water  column 

h 

m 

Depth  of  the  active  benthic  sediment  layer 

ps 

Dry  sediment  particle  density 

<l> 

fraction 

Sediment  porosity 

foe  water 

fraction 

Fraction  organic  carbon  in  solids  (water  column) 

foe  sediment 

fraction 

Fraction  organic  carbon  in  solids  (sediment) 

T 

degC 

Water  and  sediment  mean  temperature 

W 

m/sec 

Mean  wind  speed  at  10  m  above  surface 

Hyd raulic  Parameters 

X 

km 

Usage  location,  i.e.,  distance  downstream  from  the  upstream 
boundary  to  the  location  of  interest 

Geometric  option 

NA 

Select  option  of  either  entering  stream  width  and  depth  or 
entering  cross-sectional  area 

B 

m 

Stream  constant  top  width  for  option  of  entering  stream  width 
and  hydraulic  depth 

H 

m 

Stream  constant  hydraulic  depth  for  option  of  entering  stream 
width  and  hydraulic  depth 

Cross-sectional  area 
option 

NA 

For  option  of  entering  cross-sectional  area,  includes  option  of 
either  entering  a  constant  value  for  cross-sectional  area  or 
entering  a  function  related  to  flow  rate 
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Parameter  Name 

Units 

Description 

A 

m2 

Constant  cross-sectional  area  of  the  flow  for  option  of  entering 
cross-sectional  area 

a 

NA 

Parameter  in  the  function  A  =  aQh  for  option  of  entering  a 
function  for  cross-sectional  area  as  related  to  flow 

b 

NA 

Parameter  in  the  function  A  =  aQh  for  option  of  entering  a 
function  for  cross-sectional  area  as  related  to  flow 

c 

NA 

Parameter  in  the  function  H  =  cAd  for  option  of  entering  a 
function  for  hydraulic  depth  as  related  to  flow  cross-sectional 

area 

d 

NA 

Parameter  in  the  function  H  =  cAd  for  option  of  entering  a 
function  for  hydraulic  depth  as  related  to  flow  cross-sectional 

area 

Q 

m3/yr 

Constant  background  stream  flow  rate  (e.g.,  annual  mean  flow) 
at  the  head  of  the  reach  (without  any  flow  from  the  AOI) 

Constituent  Parameters 

Ci 

mg/L 

Constant  constituent  background  concentration  in  stream  at  the 
head  of  the  reach 

Cbi 

mg/kg 

Initial  constituent  concentration  in  the  sediment  bed 

kdw 

da^1 

Decay  rate  of  dissolved  phase  in  water  column 

kPw 

da^1 

Decay  rate  of  particulate  phase  in  water  column 

kdb 

da^1 

Decay  rate  of  dissolved  phase  in  sediment  bed 

kpb 

da^1 

Decay  rate  of  particulate  phase  in  sediment  bed 

Kow 

ml/ml 

Constituent  octonol-water  partition  coefficient 

Kdw 

L/kg 

Constituent  sediment-water  partition  coefficient  in  water  column 

Kdb 

L/kg 

Constituent  sediment-water  partition  coefficient  in  sediment  bed 

kv 

m/day 

Constituent  volatilization  rate 

Vd 

m/day 

Constituent  mass  transfer  rate  across  the  sediment-water 
interface  resulting  from  diffusion  of  dissolved  constituent 

MW 

g/g-mol 

Constituent  molecular  weight 

Dm 

cm2/sec 

Constituent  molecular  diffusivity  in  water 

He 

atm-m3/g-mol 

Constituent  Henry’s  law  constant 

Sedimentation  Parameters 
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Figure  16.  CMS  surface  water  model  user  interface,  first  screen  of  five. 

The  user  sets  the  model  time  step  and  the  total  simulation  time.  For  long¬ 
term  (decades  and  longer)  simulations,  the  time  step  should  be  made  as 
large  as  possible  to  reduce  run  time  without  compromising  accuracy.  For 
this  reason,  it  is  a  good  idea  to  make  runs  using  various  time  step  sizes  and 
comparing  results.  For  long-term  simulations,  o.i  to  0.2  years  is  a 
reasonable  time  step  size  to  consider.  For  short-term  simulations,  such  as  a 
month  or  less,  much  shorter  time  steps  should  be  used.  The  total  simulation 
time  is  used  to  define  how  long  of  a  period  the  model  is  run.  The  simulation 
time  should  be  long  enough  for  the  water  body  to  respond  to  the  changes  in 
loading  conditions  and  to  capture  study  needs.  For  example,  if  the  intent  is 
to  capture  conditions  during  range  use  and  years  following  it  or  future  long¬ 
term  range  use,  then  simulations  could  require  50  to  100  years. 

The  longitudinal  dispersion  coefficient  is  used  to  account  for  longitudinal 
mixing  associated  with  the  limitations  of  a  one-dimensional  assumption. 
Dispersion  accounts  for  constituent  spreading  along  the  channel  due  to 
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channel  bottom  velocity  shear  across  the  channel  and  channel  irregularities. 
Fischer  et  al.  (1979)  recommend  the  following  equation  for  estimating 
longitudinal  dispersion  Dx  (m2/sec)  in  natural  river  channels, 

TT 2  tj2 

D  =0.011 -  (68) 

Hu, 


where 

U  =  mean  flow  velocity  in  the  channel  for  the  given  flow,  m/sec 
B  =  channel  top  width  at  the  water  surface  for  the  given  flow,  m 
H  =  hydraulic  depth  of  flow,  m,  where  the  product  UHB  is  the  flow 
discharge  rate 

u*  =  shear  velocity  of  the  flow,  m/sec 

The  shear  velocity  can  be  approximated  from  0.1  U.  The  dispersion 
coefficient  should  be  estimated  for  the  annual  average  flow  rate 
(discharge)  in  the  stream. 

The  user  must  enter  the  TSS  concentration  for  the  stream.  The  version  of 
CMS  used  in  TREECS  will  maintain  a  constant  TSS  concentration  over 
time  and  throughout  the  reach.  This  version  does  not  allow  the  option  to 
conduct  suspended  solids  transport  as  the  stand-alone  version  does. 

The  active  sediment  layer  represents  the  upper  horizon  of  the  sediment 
bed  that  is  fairly  well-mixed  and  exchanges  with  the  water  column.  A 
typical  value  is  on  the  order  of  10  cm  to  20  cm. 

The  dry  sediment  density  is  typically  2650  g/L  (2.65  g/cm3).  Porosity  of  the 
active  layer  can  vary  from  about  0.5  to  0.9.  A  reasonable  value  is  0.7,  but  if 
the  sediment  layer  is  thicker  than  20  cm,  then  a  lower  value  for  porosity 
should  probably  be  used.  The  guidance  provided  for  fraction  of  organic 
carbon  within  the  RECOVERY  model  section  applies  here  too.  The  mean 
water  and  sediment  temperature  and  mean  wind  speed  should  represent 
annual  averages  if  long-term  simulations  are  being  conducted. 

Hydraulic  Parameters 

The  user  enters  the  distance  downstream  to  the  usage  (receptor)  location 
of  interest.  This  is  the  location  where  the  concentrations  must  be  known  to 
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assess  impacts  on  receptors.  Multiple  locations  can  be  assessed  in  the 
stand-alone  version  of  CMS,  but  the  TREECS  version  will  be  limited  to 
only  one  usage  location  to  reduce  system  development  complexity.  Only 
one  location  will  usually  be  needed  anyway  for  surface  water. 

As  shown  in  Figure  17,  the  user  has  the  option  of  entering  either  the  stream 
width  and  hydraulic  depth  or  the  flow  cross-sectional  area.  If  the  first  option 
is  selected,  then  the  reach-average  channel  top  width  and  hydraulic  depth 
should  be  entered  for  the  average  annual  flow  conditions.  If  the  cross- 
sectional  area  is  selected,  then  the  user  has  the  option  of  entering  a  single 
constant  area  for  the  entire  reach  or  entering  parameters  to  compute  the 
reach  area  and  depth,  given  the  flow,  which  can  change  over  time  due  to  the 
water  fluxes  (which  can  be  time-varying)  being  consumed  from  the  WFF. 
However,  for  TREECS,  only  annual  average  hydrology  is  used  presently, 
thus  the  flows  in  the  WFF  are  constant.  For  this  reason,  there  is  no  reason  to 
use  the  second  option  for  cross-sectional  area.  The  user  must  enter  the 
stream  flow  rate  Q  at  the  point  of  entry  (or  background  flow)  where  external 
loadings  (such  as  runoff  from  AOI)  are  entering,  which  is  the  beginning  of 
the  modeled  stream  reach.  As  stated  above,  this  background  flow  does  not 
vary  in  time. 


Figure  17.  The  Hydraulic  Parameters  user  input  screen  of  CMS. 
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Constituent  Parameters 

Constituent  parameters  must  be  entered  for  each  modeled  constituent  on 
the  constituent  parameters  screen,  which  is  not  shown  here.  The  user  must 
enter  the  background  stream  concentration,  which  is  the  concentration  in 
the  stream  at  the  head  of  the  reach  before  mixing  with  the  AOI  loading. 
Typically,  this  value  can  be  set  to  zero  unless  there  is  a  known  background 
concentration.  The  initial  concentration  of  MC  in  the  sediment  bed  must  be 
entered.  However,  it  may  be  set  to  zero  in  many  applications  since  this  input 
can  have  little  significance  for  long-term  simulations  with  continuing 
loadings  since  the  initial  bed  concentration  will  eventually  change  over  time 
to  reflect  the  loadings.  However,  if  the  user  is  not  conducting  a  long-term 
analysis  with  continuing  loadings  and  needs  to  input  initial  concentrations 
in  sediment  and  water,  he  or  she  may  do  so.  There  is  no  input  for  initial 
water  concentration  since  streams  are  usually  flushed  rather  quickly  of  any 
initial  conditions  for  the  water  column. 

The  four  decay  rates  can  be  set  to  zero  for  conservative  predictions  of 
constituent  fate.  Alternatively,  the  user  can  input  rates  for  dissolved  and 
particulate  constituent  phases  for  both  the  water  column  and  the  sediment 
bed.  Little  guidance  can  be  provided  for  values  to  use  since  rates  can  vary 
widely  among  constituent  chemicals,  as  well  between  water  and  sediment 
and  particulate  and  dissolved.  The  assumption  is  often  made  that  the 
particulate  phase  constituent  does  not  decay,  while  dissolved  phase  does. 
Some  chemicals  may  decay  faster  in  the  water  column  than  in  the  bed, 
while  others  behave  the  opposite.  It  is  noted  that  decay,  or  degradation,  in 
this  model  includes  all  loss  types  (biodegradation,  photolysis,  hydrolysis, 
and  oxidation)  except  volatilization,  which  is  a  separate  process  as  in  the 
RECOVERY  model.  The  decay  rate  for  metals  is  usually  set  to  zero. 

The  user  has  the  option  of  choosing  either  to  input  the  Octonol -water 
partition  coefficient,  which  is  used  by  the  user  interface  to  compute  Kd  for 
organic  chemicals,  or  input  the  Kd  for  the  water  column  and  sediment  bed. 
The  second  option  should  be  chosen  for  non-organic  constituents,  such  as 
metals,  but  can  also  be  used  for  any  constituent  if  the  user  has  better 
information  on  Kd.  If  a  value  for  Kow  is  in  the  constituent  database  and  that 
option  is  selected,  that  value  is  automatically  passed  to  the  input  field  for 
this  parameter.  The  model  computes  Kd  when  Kow  is  provided  using  the 
same  methods  as  described  for  the  RECOVERY  model. 
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The  user  can  either  have  the  user  interface  compute  the  rate  parameters 
for  volatilization  and  dissolved  phase  mass  transfer  between  water  column 
and  bed,  or  the  user  can  input  these  two  rates.  The  methods  for  these 
calculations  are  described  by  Fant  and  Dortch  (2007),  and  all  of  the 
information  required  for  the  calculations  is  contained  within  the  interface 
input. 

Values  for  the  three  chemical-specific  properties,  molecular  weight, 
molecular  diffusivity  in  water,  and  Henry’s  law  constant,  are  automatically 
passed  from  the  constituent  database  if  values  are  present.  Otherwise,  the 
user  must  enter  the  missing  values. 

Sedimentation  Parameters 

Similar  to  the  RECOVERY  model,  the  version  of  CMS  in  TREECS  requires 
the  user  to  enter  two  of  the  three  sedimentation  parameters,  settling  rate, 
burial  rate,  and  resuspension  rate,  and  the  remaining  unspecified  parameter 
is  computed  from  a  steady-state  solids  balance.  The  three  rates  are  constant 
throughout  the  study  reach  and  over  time.  It  is  usually  easier  to  estimate  the 
settling  and  resuspension  rates  for  rivers.  The  user  must  ensure  that  a 
positive  burial  rate  is  computed.  Estimation  of  settling  and  resuspension 
rates  is  beyond  the  scope  of  this  report,  but  there  is  a  large  body  of 
knowledge  regarding  the  estimation  of  these  parameters.  Settling  rate  can 
be  estimated  from  TSS  mean  particle  size  and  Stokes  law.  Estimation  of 
resuspension  rates  is  discussed  by  Fant  and  Dortch  (2007). 

Output 

The  output  of  the  CMS  includes  the  following  for  each  MC,  which  are 
viewable  via  a  special  viewer  developed  for  CMS: 

•  Water  column  total  or  dissolved  concentration  (mg/L)  versus  stream 
length  for  given  times 

•  Time  animation  of  water  column  total  or  dissolved  concentration 
(mg/L)  versus  stream  length 

•  Bed  sediment  total  concentration  on  a  sediment  mass  basis  (mg/kg) 
versus  stream  length  for  given  times 

•  Time  animation  of  bed  sediment  total  concentration  on  a  sediment 
mass  basis  (mg/kg)  versus  stream  length 

•  Water  column  total  or  dissolved  concentration  (mg/L)  versus  time  for 
given  distances  along  the  stream 
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•  Distance  animation  of  water  column  total  or  dissolved  concentration 
(mg/L)  versus  time 

•  Bed  sediment  total  concentration  on  a  sediment  mass  basis  (mg/kg) 
versus  time  for  given  distances  along  the  stream 

•  Distance  animation  of  bed  sediment  total  concentration  on  a  sediment 
mass  basis  (mg/kg)  versus  time 

In  addition  to  the  CMS  viewer,  there  are  three  general  types  of  FRAMES 

viewers  that  will  provide  output  results.  These  include: 

•  SCF  (Sediment  Concentration  File)  graphical  view  -  time  series  of  bed 
sediment  total  concentration  on  a  sediment  mass  basis  (mg/kg)  and 
bed  sediment  dissolved  concentration  on  a  water  volume  basis  (pore 
water  concentration,  mg/L) 

•  WCF  graphical  view  -  time  series  of  water  column  total  concentration 
(mg/L)  and  dissolved  concentration  (mg/L) 

•  WFF  graphical  view  -  time  series  of  water  column  particulate 
(adsorbed)  mass  flux  (g/yr)  and  dissolved  mass  flux  (g/yr)  associated 
with  the  constituents  in  the  steam  flow 
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8  Other  Considerations 

There  are  several  additional  features  that  merit  discussion.  These  features 
include:  soil  interflow  and  groundwater  discharge  to  surface  water;  fate  of 
mixtures;  fate  of  degradation  products;  sensitivity  and  uncertainty,  and 
modeling  fate  of  water  miscible  constituents.  Each  of  these  features  is 
discussed  below  as  well  as  plans  for  how  each  will  be  handled  in  TREECS 
Tier  2. 

Soil  Interflow  and  Groundwater  Discharge 

Interflow 

Interflow  through  the  vadose  zone  or  soil  to  surface  water  is  a  minor 
pathway  in  most  cases,  but  it  could  be  a  potential  pathway  for  surface  soils 
having  a  high  hydraulic  conductivity  with  an  impermeable  or  semi- 
impermeable  soil  layer  at  a  shallow  soil  depth,  thus  creating  a  perched  water 
table.  Since  the  TREECS  Tier  1  and  2  models  are  based  on  average  annual 
hydrology,  interflow  must  be  added  through  specification  of  the  fraction  of 
average  annual  infiltration  that  is  lost  to  interflow  before  entering  the 
vadose  zone  model.  Interflow  is  caused  by  infiltration  flow  that  is  greater 
than  the  maximum  percolation  rate  (i.e.,  groundwater  recharge  rate)  of  the 
vadose  zone  layer,  which  is  the  saturated  hydraulic  conductivity,  Ks.  A 
conceptual  schematic  of  flow  within  the  soil,  vadose,  aquifer  system  is 
shown  in  Figure  18.  After  prompting  the  user  to  enter  a  value  of  Ks  for  the 
vadose  zone  layer,  the  soil  model  user  interface  will  automatically  provide 
an  estimate  of  the  interflow  fraction  as  follows.  If  the  soil  infiltration  rate, 
qw,  is  less  than  or  equal  to  Ks,  the  fraction  of  interflow,  Fif,  is  estimated  to  be 

d  _ 

zero.  If  qw  is  greater  than  Ks,  then  Fif  =  — - - .  The  user  can  also  inter 

Qw 

his/her  own  estimate  for  Fif.  In  many  cases,  the  interflow  fraction  should  be 
set  to  zero,  which  will  send  all  of  the  infiltration  water  to  the  vadose  zone 
and  aquifer  models. 

No  fate  processes,  such  as  sorption  or  degradation,  are  applied  to  the 
interflow  flux.  Thus,  the  fraction  of  interflow  will  be  multiplied  times  the 
infiltrating  flow  rate  and  mass  fluxes  computed  by  the  soil  model  to 
produce  the  interflow  water  flow  rate  and  mass  fluxes.  These  values  will 
then  be  added  to  the  surface  runoff  flow  rate  and  dissolved  mass  fluxes 
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Figure  18.  Conceptual  schematic  of  surface  and  sub-surface 
hydrology. 


computed  by  the  soil  model,  and  the  combined  results  will  be  written  to 
the  WFF  surface  water.  The  soil  model  must  also  reduce  the  water  flow 
and  mass  fluxes  to  vadose  zone  by  the  same  amount  before  writing  to  the 
WFF  vadose  file. 

Groundwater  Discharge 

During  dry  seasons  when  surface  stream  flows  are  low,  groundwater 
discharge  can  contribute  a  major  portion  of  the  stream  flow.  For  some 
watershed  conditions,  groundwater  discharge  can  be  the  primary 
contributor  to  flow  in  small  streams.  For  these  reasons,  it  is  necessary  to 
include  a  contribution  of  MC  loading  from  groundwater  flow  in  surface 
water. 

A  module  must  be  developed  to  provide  groundwater  discharge  to  surface 
water.  The  new  module  will  consume  the  MEPAS  Aquifer  water  flux  file 
(WFF)  which  contains  aquifer  flow  (m3/yr)  and  aquifer  mass  flux  (g/yr)  at  a 
flux  location,  which  is  specified  in  the  MEPAS  Aquifer  model  user  interface. 
The  flux  location  in  this  case  should  be  the  distance  from  the  center  of  the 
AOI  to  the  point  in  the  landscape  where  the  stream  or  water  body  and 
aquifer  have  crossed  in  the  landscape.  If  the  target  surface  water  point  of 
interest  (i.e.,  receptor  target  location)  is  further  downstream,  then  there  are 
no  further  mass  losses  or  mass  flux  attenuation  as  water  travels  from  the 
point  of  groundwater  intersection  to  the  water  body  target  location. 
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Since  the  receiving  surface  water  body  models  within  TREECS  are 
designed  to  receive  surface  water  runoff  and  fluxes  from  the  AOI  soil,  the 
new  module  must  also  be  able  to  combine  surface  and  aquifer  water  flux 
files  (WFFs)  into  a  combined  surface  water  WFF.  A  Plus  operator  was 
developed  for  ARAMS  that  allows  one  to  combine  multiple  WFFs  into  one, 
but  the  WFFs  had  to  be  of  the  same  type,  such  as  aquifer  or  surface  water. 
The  reason  for  this  is  that  the  surface  water  WFF  contains  mass  flux  for 
dissolved  and  adsorbed  mass  flux,  whereas  the  aquifer  WFF  contains  only 
dissolved  mass  flux. 

The  new  module  must  serve  two  functions.  It  must  perform  an  operation 
on  the  aquifer  WFF  for  groundwater  discharge  to  surface  water,  and  it 
must  allow  combining  of  surface  water  and  aquifer  WFFs.  The  new 
module  is  referred  to  here  as  the  Plus  operator  for  Surface  water  and 
Groundwater,  or  Plus-SG  for  short. 

The  Plus-SG  module  should  have  a  brief  user  interface  that  will  request  the 
user  to  input  the  groundwater  flow  rate  entering  the  receiving  surface 
water  body.  The  flow  rate  input  units  can  vary,  but  the  flow  rate  will  be 
converted  to  m'Vyr  to  be  consistent  with  WFF  specifications.  The 
groundwater  discharge  to  surface  water  will  be  assumed  to  be  constant 
over  time,  which  is  fairly  consistent  with  base  flow  conditions  in  small 
streams.  The  Plus-SG  module  will  divide  the  aquifer  WFF  mass  flux  by  the 
aquifer  WFF  flow  at  the  flux  location  to  obtain  an  aquifer  concentration  at 
the  location  of  groundwater  discharge  to  surface  water.  This  concentration 
will  be  multiplied  by  the  user-input  groundwater  discharge  rate  to  obtain 
values  for  mass  flux  (g/yr)  from  groundwater  to  surface  water.  The  water 
flux  rate  will  be  the  constant  value  of  groundwater  discharge  to  surface 
water  input  by  the  user. 

The  groundwater  discharge  water  and  mass  fluxes  must  be  combined  with 
the  surface  water  WFF.  The  surface  water  WFF  dissolved  mass  fluxes  that 
enter  the  Plus-SG  will  be  combined  with  the  groundwater-to-surface-water 
mass  flux  to  obtain  a  total  dissolved  mass  flux  entering  the  water  body. 

The  surface  water  WFF  particulate  (adsorbed)  mass  fluxes  that  enter  the 
Plus-SG  will  be  passed  unaltered  to  the  Plus-SG  output  WFF,  which  will  be 
consumed  by  the  water  body.  The  surface  water  WFF  water  flux  will  be 
combined  with  the  groundwater-to-surface-water  flow  to  obtain  the  total 
flow  entering  the  water  body.  Thus,  the  Plus-SG  operator  will  output  a 
WFF  that  is  similar  to  the  surface  water  WFF,  except  that  surface  water 
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and  groundwater-to-surface-water  and  mass  fluxes  have  been  combined. 
This  output  file  will  be  consumed  by  the  surface  water  models. 

The  existing  Plus  Operators  are  able  to  handle  time  series  with  different 
time  points  when  combining  fluxes  and  flows.  This  feature  must  be 
preserved  in  the  Plus-SG  module. 

The  RECOVERY  model  uses  only  the  WFF  mass  flux  or  load  (g/yr)  as 
input;  water  flux  from  WFF  is  not  used.  Since  the  RECOVERY  model 
assumes  a  fully  mixed  water  column,  it  does  not  matter  that  the  mass 
fluxes  for  surface  water  and  groundwater  are  combined.  The  user  must 
enter  a  water  body  flow  through  rate  in  the  model  user  interface  since 
other  water  sources  could  enter.  This  flow  rate  should  represent  the 
combined  flows  from  AOI  runoff,  groundwater  discharge,  and  other 
sources  of  runoff  not  associated  with  the  firing  range  or  AOI. 

The  CMS  uses  both  the  WFF  load  and  water  flux  (flow  rate)  as  input.  The 
CMS  does  not  presently  have  a  feature  to  allow  for  flow  distribution  along 
the  river  reach,  as  might  occur  with  groundwater  discharge.  Thus,  all  of  the 
flow  from  AOI  surface  water  runoff  and  groundwater  enter  at  the  head  of 
the  modeled  reach.  The  user  can  also  input  a  constant  background  flow, 
which  is  flow  entering  the  head  of  the  reach  before  the  flows  from  AOI 
runoff  and  groundwater  discharge  are  added.  Consideration  should  be  given 
to  adding  a  distributed  flow  feature  in  the  future  for  groundwater  discharge. 

Fate  of  Mixtures 

Some  MC  can  exist  in  mixtures,  such  as  several  HE,  including  Composition 
(Comp)  B,  Octol,  and  Tritonal,  which  are  mixtures  of  RDX  and  TNT,  HMX 
and  TNT,  and  TNT  and  aluminum,  respectively.  Mixtures  can  affect  the 
dissolution  of  individual  constituents  within  the  mixture. 

Taylor  et  al.  (2009a)  found  that  for  Comp  B,  the  dissolution  rate  of  less 
soluble  RDX  controlled  the  dissolution  rate  of  more  soluble  TNT.  The  TNT 
dissolution  rate  was  essentially  the  RDX  dissolution  rate  times  the  initial 
bulk  RDX/TNT  mass  ratio.  However,  initially,  the  dissolution  rates  of  RDX 
and  TNT  were  about  the  same  without  the  mass  ratio  multiplier.  As 
discussed  in  Appendix  B,  the  TREECS  dissolution  model  accurately 
predicted  Comp  B,  RDX,  and  TNT  dissolution  masses  using  the  solubility  of 
the  mixture  and  the  respective  mass  fractions  of  the  mixture,  suggesting 
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that  the  solubility  of  the  mixture  could  be  the  proper  approach  for  applying 
the  TREECS  dissolution  model. 

Taylor  et  al.  (2009a  and  2009b)  found  that  the  dissolution  rate  of  the  TNT 
in  Tritonal  occurred  at  the  pure  TNT  dissolution  rate  as  if  the  aluminum 
flakes  had  no  influence.  However,  the  AF  fitting  parameter  in  their  model 
had  to  be  reduced  from  2  for  pure  TNT  to  1.22  for  Tritonal,  or  nearly  by  half 
as  discussed  in  Appendix  B.  If  the  solubility  of  TNT  is  reduced  by  half,  the 
TREECS  dissolution  model  also  matches  the  observed  TNT  dissolution 
results  fairly  well,  which  is  also  discussed  in  Appendix  B.  Taylor  et  al. 
(2009a)  found  that  the  dissolution  rate  of  less  soluble  HMX  did  not  fully 
restrain  the  dissolution  rate  of  TNT  for  Octol,  where  TNT  dissolution 
occurred  at  a  rate  between  the  independent  TNT  and  HMX  controlled  rates. 

The  above  discussion  suggests  that  the  mechanisms  for  dissolution  of 
mixtures  is  still  not  well  understood  and  merits  further  study.  Additionally, 
MC  residue  can  occur  as  a  result  of  munitions  that  use  both  pure  explosive 
components  and  mixtures,  such  as  pure  RDX  and  Comp  B.  The  MIDAS 
database  does  not  provide  information  on  munitions  mixtures,  such  as  the 
amount  of  Comp  B;  rather,  it  only  provides  the  mass  of  each  explosive 
component,  such  as  RDX  and  TNT.  As  a  result  the  amount  of  mixture 
residue  may  not  be  known  or  could  be  difficult  to  estimate.  To  account  for 
mixtures,  it  is  also  necessary  to  conduct  a  mass  balance  on  each  mixture 
residue  in  addition  to  the  pure  constituent  mass  balances.  These  limitations 
and  requirements  make  it  difficult  to  model  the  fate  of  mixtures.  Inclusion 
of  mixtures  at  this  time  could  result  in  delays  in  delivering  an  initial 
TREECS  Tier  2  capability.  Therefore,  it  is  recommended  that  explicit 
treatment  of  mixtures  not  be  included  in  Tier  2  at  this  time,  but  this 
capability  should  be  a  consideration  for  future  versions. 

Accordingly,  all  MC  will  be  assumed  to  be  in  pure  form,  and  the 
dissolution  rates  will  be  computed  with  Equations  40-43  using  pure 
constituent  properties.  This  may  not  be  a  bad  assumption  for  the  less 
soluble  constituents  if  their  dissolution  rate  really  controls  the  dissolution 
rate  for  mixtures.  This  is  a  conservative  assumption  for  the  more  soluble 
constituents  since  they  will  be  predicted  to  dissolve  faster  than  could 
actually  occur  for  mixtures,  thus  promoting  their  mobility  in  water 
systems.  Should  mixtures  be  included  in  future  versions  of  TREECS,  it  is 
recommended  at  this  time  that  the  solubility  of  the  mixture  be  used  for  the 
dissolution  modeling  if  that  solubility  can  be  estimated. 
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Fate  of  Degradation  Products 

The  degradation  (or  daughter)  products  of  some  MC,  such  as  TNT,  can 
pose  more  concern  than  the  parent  compound.  The  fate  of  each  daughter 
product  of  concern  should  be  tracked  through  each  media.  This  requires 
modeling  the  fate  and  transport  of  both  parent  MC  and  daughter  products 
within  each  medium  model  where  the  source  of  the  daughter  product  is 
the  degradation  flux  of  the  parent  constituent.  Presently,  the  TREECS 
system  and  several  of  the  models  (soil  and  surface  water)  do  not  provide 
the  means  to  model  degradation  products.  Providing  such  a  capability  will 
require  extensive  development  and  modifications.  Thus,  the  initial 
versions  of  TREECS  will  not  have  this  capability,  but  this  extension  is 
recommended  for  future  versions. 

Sensitivity  and  Uncertainty 

The  capability  to  assess  parameter  and  input  sensitivity  and  output 
uncertainty  (S/U)  will  be  included  in  Tier  2,  but  not  Tier  1,  of  TREECS.  This 
capability  will  utilize  modules  and  specifications  that  exist  within 
ARAMS™.  The  unaltered  FRAMES  S/U  module  within  ARAMS™  will  be 
used  in  Tier  2  for  setting  up  inputs  for  the  assessment.  The  methods  are 
based  on  Monte  Carlo  simulation  using  Latin  Hypercube  sampling  of 
random  inputs.  The  user  specifies  the  type  of  distribution  (such  as  normal, 
log  normal,  etc.)  and  the  distributions  statistics  for  each  declared  uncertain 
input.  The  output  can  be  sampled  to  provide  the  occurrence  frequency  of 
specific  measures,  such  as  peak  concentration,  or  to  provide  uncertainty 
confidence  bands  on  output  measures,  such  as  concentration  time  series. 

Modeling  Fate  of  Water  Miscible  Constituents 

Some  MC  may  be  water  miscible  or  easily  and  quickly  dissolved  in  water, 
such  as  perchlorate  constituents.  There  was  a  question  of  whether  the  Tier 
2  soil  model  would  be  able  to  handle  the  very  rapid  dissolution  that  occurs 
with  such  highly  soluble  MC.  In  such  cases,  the  dissolution  flux  Fctis  in 
Equation  2  rapidly  approaches  the  MC  loading  rate  L,  the  solid  phase 
erosion  flux  Fes  rapidly  approaches  zero,  and  the  solid  phase  mass  Ms  is 
small.  Thus,  water  miscible  MC  mostly  exists  as  non-solid  phase  mass  Mns- 
Several  constituent  properties  are  no  longer  important  for  miscible  MC, 
including  initial  solid  phase  concentration,  initial  particle  diameter  dsm 
and  solid  phase  particle  mass  density  psm.  Tests  were  run  to  determine  the 
response  of  the  Tier  2  soil  model  for  miscible  MC. 
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A  spreadsheet  version  of  the  Tier  2  soil  model  was  applied  for  a  case  of 
water  miscible  MC  using  potassium  perchlorate  as  the  MC,  which  has  a 
solubility  of  about  20,000  mg/L.  The  initial  particle  diameter  was  set  to 
500  pm.  The  spreadsheet  model  used  a  first-order  accurate  Euler  time 
integration  scheme  to  solve  the  differential  equations.  With  a  time  step  of 
0.01  years,  the  model  did  not  converge  on  a  correct  solution.  With  a  time 
of  0.001  years,  the  model  did  properly  converge  to  values  fairly  close  to  the 
correct  answers.  For  smaller  initial  particle  diameters,  dissolution  was 
faster,  so  smaller  times  steps  were  required  to  keep  a  stable  solution. 

The  Tier  2  soil  model  was  also  programmed  using  the  Visual  Basic 
language  with  two  numerical  solution  options.  One  option  is  the  regular 
fourth-order  Runge-Kutta  integration  method,  which  uses  a  constant  time 
step.  The  other  option  is  the  adaptive  Runge-Kutta-Fehlberg  (Press  et  al. 
1996)  integration  method,  which  allows  the  time  step  to  adjust  as 
necessary  to  maintain  a  stable  solution.  The  programmed  model  was  run 
for  the  same  test  case  of  perchlorate  as  was  run  for  the  spreadsheet  version 
with  an  initial  particle  diameter  of  500  pm.  With  a  constant  time  step  of 
0.01  years,  the  model  ran  to  a  stable,  but  incorrect  solution.  With  a 
constant  time  step  of  0.001  years,  the  model  ran  to  a  stable  and  correct 
solution  with  nearly  the  same  results  as  the  spreadsheet  model.  The  model 
was  also  run  with  the  auto-time-stepping  option.  That  run  converged  to 
the  correct  solution  without  having  to  know  the  time  step  required  for  a 
valid  solution.  However,  when  the  model  was  run  with  an  initial  particle 
diameter  of  1.0  pm  and  with  auto-time-stepping,  the  model  had  much 
difficulty  finding  a  small  enough  time  step  to  complete  the  run. 

Given  the  potential  stability  and  accuracy  problems  noted  above,  it  was 
decided  that  a  miscible  MC  feature  must  be  added  to  the  Tier  2  soil  model. 
An  additional  input  variable  is  needed  to  declare  if  an  MC  is  miscible.  If  an 
MC  is  miscible,  then  the  dissolution  flux  Fdis  is  set  equal  to  the  loading  flux  L 
for  each  time  step  during  the  model  solution.  Additionally,  the  calculations 
for  particle  diameter  dsm  and  particle  specific  surface  area  a  are  not  required 
for  miscible  MC,  and  the  solid  phase  erosion  flux  Fes  should  be  zero.  Inputs 
for  initial  particle  diameter,  initial  solid  phase  concentration,  and  particle 
density  are  not  required  for  miscible  MC.  Testing  of  the  miscible  MC  feature 
showed  that  the  above  changes  worked  very  satisfactorily  without  any 
stability  or  accuracy  issues. 
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9  Summary 

TREECS  Tier  2  assessment  methods  will  require  more  detailed  site  data, 
and  will  require  more  knowledge  and  skill  to  apply  than  Tier  l,  but  Army 
environmental  staff  that  has  a  cursory  understanding  of  multi-media  fate 
and  transport  should  have  no  problem  applying  Tier  2.  The  Tier  2  approach 
will  allow  time-varying  analyses  with  solid  phase  dissolution  and  system 
losses,  such  as  degradation  and  volatilization.  Such  analyses  should  provide 
more  accurate  predictions  with  generally  lower  receiving  media 
concentrations  due  to  mediating  effects  of  transport  phasing/dampening 
and  natural  attenuation/degradation.  Tiers  l  and  2  focus  on  contaminant 
stressors  and  human  and  ecological  health  end  point  metrics.  This  report 
describes  the  approach  and  formulations  that  will  be  used  for  Tier  2  of 
TREECS. 

Like  Tier  l,  Tier  2  addresses  impact  areas  of  interest  and  does  not  consider 
firing  points  at  this  time.  Also  like  Tier  l,  Tier  2  will  compare  forecasted 
media  concentrations  against  protective  human  and  ecological  health 
benchmark  concentrations.  Receiving  media  include  groundwater,  surface 
water,  and  surface  water  sediments.  Additionally,  Tier  2  uses  average 
annual  hydrologic  input  like  Tier  l.  However,  Tier  2  should  be  more 
accurate  than  Tier  l  for  assessing  MC  fate,  especially  for  metals,  since  the 
time  domain  and  dissolution  are  included.  Tier  2  should  be  used  when 
Tier  l  indicates  that  a  protective  health  benchmark  will  be  exceeded. 

The  Tier  2  approach  for  TREECS  is  based  on  solving  mass  balance 
equations  for  MC  mass  that  evolves  over  time.  The  soil  model  assumes  that 
a  single  homogeneous  layer  or  soil  compartment  represents  the  AOI  with 
spatially  uniform  properties  and  MC  concentrations.  However,  MC  residue 
loadings  to  soil  and  computed  soil  MC  mass,  concentrations  and  fluxes 
(losses,  gains,  and  exports)  can  vary  over  time.  Including  time  as  a 
dimension,  which  was  not  included  for  Tier  l,  allows  consideration  of  both 
solid  and  non-solid  phase  MC  mass  with  dissolution  of  the  solid  phase 
transferring  mass  to  the  non-solid  phase.  The  non-solid  phase  mass  consists 
of  aqueous  dissolved,  aqueous  adsorbed-to-soil,  and  gas  or  vapor  in  air 
components.  Equilibrium  portioning  among  the  non-solid  components  is 
assumed,  which  collapses  the  non-solid  mass  balance  into  a  single 
differential  equation.  Thus,  the  Tier  2  soil  model  solves  two  ordinary 
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differential  equations,  one  for  solid  phase  and  one  for  the  non-solid  phase 
mass  balance. 

Tier  2  receiving  media  models  consist  of  the  MEPAS  vadose  zone  model, 
MEPAS  aquifer  model,  RECOVERY  surface  water  model,  and  CMS  surface 
water  model.  Both  surface  water  models  include  benthic  sediments.  MC 
mass  can  move  from  soil  to  vadose  zone,  and  from  vadose  zone  to  aquifer. 
MC  mass  can  also  move  from  soil  to  surface  water.  The  user  decides  which 
surface  water  model  to  use,  either  RECOVERY  or  CMS.  MC  mass  can  move 
as  interflow  from  the  soil  compartment  to  surface  water  where  interflow  is 
the  excess  infiltration  that  cannot  percolate  through  the  vadose  zone.  There 
are  no  mass  losses  or  gains  associated  with  the  interflow.  Additionally,  MC 
mass  can  move  from  groundwater  to  surface  water.  The  initial  version  of 
Tier  2  will  not  allow  MC  mass  to  move  from  surface  water  to  surface  water, 
such  as  transport  from  a  stream  to  a  lake,  although  this  feature  could  be 
added  later  if  necessary.  The  RECOVERY  model  is  better  suited  for 
modeling  standing  or  pooled  water,  such  as  lakes,  ponds,  wetlands,  etc, 
whereas  the  CMS  is  better  suited  for  modeling  flow  in  streams  and  rivers, 
including  tidal  rivers. 

The  MC  residue  mass  loading  model  is  described.  The  residue  mass  loading 
model  takes  munitions  use  information  and  develops  mass  loadings  for 
each  MC.  The  MIDAS  database  is  used  to  establish  the  amount  of  MC  within 
each  munitions  type.  Sources  of  MC  residue  within  the  loading  model 
include  unexploded  mass  for  high-order  and  low-order  detonations  and 
sympathetic  detonation  of  duds.  Corrosion  and  release  of  MC  inside  duds  is 
not  considered  at  this  time.  The  number  of  rounds  fired  per  year  must  be 
entered  for  each  munitions  item  used.  Estimates  of  yield  percentages  for 
low-order,  high-order,  and  sympathetic  dud  detonations  and  occurrence 
percentages  of  duds,  low-order  detonations,  and  sympathetic  dud 
detonations  must  be  entered. 

The  fate  formulations  within  the  Tier  2  soil  model  are  developed  and 
presented,  and  each  of  the  receiving  media  fate  models  is  described  as 
well.  The  input  requirements  for  each  fate  model  are  presented  and 
information  for  obtaining  inputs  is  provided.  The  output  from  each  fate 
model  is  also  defined. 

Some  in-depth  study  and  analysis  were  required  for  establishing 
appropriate  mathematical  descriptions  for  several  fate  processes,  including 
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dissolution,  volatilization,  and  solid  phase  particle  erosion.  Each  of  the 
processes  can  be  quite  complex.  A  simplified  volatilization  formulation 
provides  reasonable  flux  estimates  without  requiring  a  more  complex  model 
of  depth-varying  soil  properties  and  concentrations.  A  generic  dissolution 
formulation  is  developed  that  is  relatively  easy  to  apply  and  does  not  require 
calibration  parameters.  The  dissolution  formulation  agrees  with  previous 
formulations  and  field  measurements.  Bed  load  formulations  were  explored 
for  describing  solid  phase  MC  erosion.  The  Einstein  and  Brown  formulation 
seemed  to  perform  favorably  for  estimating  soil  erosion,  and  it  can  probably 
be  applied  for  mixtures  of  soil  and  solid  phase  MC.  However,  a  much 
simpler  approach  is  recommended  for  Tier  2  which  is  based  on  the  USLE. 
The  Einstein  and  Brown  approach  may  be  more  appropriate  for  short-term, 
storm-event  modeling  of  MC  fate. 

A  Plus-SG  Operator  module  must  be  developed  to  allow  import  fluxes  of 
both  surface  water  erosion/runoff  and  groundwater  discharge  into 
receiving  surface  water.  This  operator  will  combine  the  two  types  of  fluxes 
after  requiring  that  the  user  input  the  groundwater  discharge  flow  rate. 

The  fate  of  mixtures  is  affected  by  the  dissolution  rate  of  mixtures,  which 
still  remains  a  topic  for  further  research.  Proper  modeling  of  mixtures 
requires  conducting  mass  balances  on  each  mixture  mass,  as  well  as  mass 
balances  for  each  constituent  in  the  mixture.  The  mixture  mass  residue 
loading  rates  onto  the  range  must  also  be  known  or  estimated.  This 
requires  knowing  which  mixtures  are  contained  in  munitions.  The  MIDAS 
database  does  not  provide  this  information;  rather,  it  provides  the  mass  of 
each  constituent  or  explosive  pure  component.  Therefore,  dissolution  of 
mixtures  will  not  be  implemented  into  Tier  2  at  this  time  although  this 
feature  should  be  considered  for  future  versions  of  TREECS. 

The  capability  to  simulate  the  fate  of  degradation  products  also  will  not  be 
implemented  into  Tier  2  at  this  time.  This  feature  will  require  additional 
development  for  the  soil  model  as  well  as  extensive  modifications  to  the 
TREECS  system  and  the  receiving  water  models.  Such  development  would 
preclude  the  release  of  Tier  2  during  the  funding  cycle  of  the  present 
TREECS  project.  This  capability  should  be  considered  for  future  versions 
of  TREECS. 

Tier  2  will  consider  only  one  initial  particle  size  for  solid  phase  MC.  The 
Tier  2  soil  model  will  compute  the  time  evolution  of  the  single  particle  size, 
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which  affects  the  time-varying  dissolution  rate.  Particle  size  also  affects 
the  solid  phase  particle  and  soil  mixture  erosion  rate  if  a  bed  load  type 
function  is  used  for  erosion.  The  inclusion  of  multiple  particle  sizes 
requires  conducting  a  mass  balance  on  each  particle  size  class  for  each  MC. 
This  feature  would  also  require  input  of  the  distribution  of  initial  solid 
phase  particle  sizes  for  each  MC.  Including  multiple  solid  phase  particle 
sizes  substantially  increases  the  input  data  requirements  and  modeling 
requirements;  thus,  this  feature  should  be  considered  carefully  before 
deciding  to  implement  it  for  future  versions  of  TREECS. 

The  capability  to  assess  parameter  and  input  sensitivity  and  output 
uncertainty  (S/U)  will  be  included  in  Tier  2,  but  not  Tier  1,  of  TREECS. 
This  capability  will  utilize  modules  and  specifications  that  exist  within 
ARAMS™.  The  output  can  be  sampled  to  provide  the  occurrence 
frequency  of  specific  measures,  such  as  peak  concentration,  or  to  provide 
uncertainty  confidence  bands  on  output  measures,  such  as  concentration 
time  series.  The  S/U  module  uses  Monte  Carlo  simulation  with  Latin 
Hypercube  sampling. 

The  Tier  2  soil  model  will  handle  miscible  MC,  such  as  perchlorate  that 
dissolves  rapidly  in  water,  without  any  special  consideration.  However,  the 
adaptive  time  step  integration  scheme  should  be  the  default  method  for 
the  final  programmed  Tier  2  soil  model  in  TREECS. 

The  TREECS  advisory  panel  expressed  an  interest  in  the  development  of 
an  additional  TREECS  tool  for  assessing  single  storm  events  and  first  flush 
MC  concentrations  in  streams.  This  feature  was  not  included  in  this  initial 
version  of  Tier  2  in  TREECS.  However,  it  is  feasible  to  add  this  feature 
sometime  in  the  future. 
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Appendix  A:  Solid  Phase  MC  Particle  Erosion 

It  is  feasible  for  solid  phase  MC  particles  to  erode  with  soil  during  rainfall 
events,  washing  into  streams  and  traveling  to  target  water  bodies.  In  fact, 
results  of  laboratory  studies  (Larson  et  al.  2004  and  2005)  indicate  that 
particulate  lead  in  runoff  is  a  major  pathway  for  lead  migration  from  firing 
range  soils.  However,  these  studies  did  not  determine  whether  the 
particulate  lead  was  soil  adsorbed  aqueous  phase  or  solid  phase  lead 
particles.  Additionally,  laboratory  studies  of  erosion  may  not  reflect  the 
true  nature  of  what  happens  in  the  field  where  overland  flow  and  eroded 
material  must  travel  much  greater  distances. 

Soil  adsorbed  aqueous  phase  transport  is  already  taken  into  account  in  the 
Tier  2  soil  model  as  described  in  Chapter  4.  Solid  phase  MC  erosion  was 
identified  as  a  potential  pathway  in  Chapter  4,  but  methods  for  computing  it 
were  not  defined.  Various  sediment  transport  formulations  were  reviewed 
in  an  attempt  to  find  a  way  to  describe  the  process  of  solid  phase  MC 
particle  erosion.  This  process  is  quite  complex  since  the  size  and  specific 
gravity  of  soil  and  MC  particles  can  affect  the  erosion  rate.  Thus,  it  is 
possible  that  the  properties  of  the  mixture  of  soil  and  MC  determine  erosion 
rates,  depending  on  the  quantity  of  MC  in  the  soil. 

Julien  (1995)  presents  various  methods  for  estimating  bed  load  transport 
in  channels.  The  methods  were  reviewed,  and  two  methods  (Duboy‘s 
equation  and  the  Einstein  and  Brown  equations)  were  evaluated  to  explore 
their  potential  use.  The  two  methods  were  applied  to  hypothetical 
conditions  as  well  as  laboratory  data  presented  by  Larson  et  al.  (2004). 
Although  the  details  of  the  two  methods  and  their  application  results  are 
not  presented  here,  the  general  conclusions  drawn  are  described  below. 

Duboy’s  Equation 

The  Duboy’s  equation  is  stated  as, 

0.173  ,  v 

Qbv  70.75  To  \To  —T-c) 

(Al) 


where 
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qbv  =  volume  of  bed  material  load  per  unit  width,  ft2/sec 
ds  =  particle  diameter,  ft 
To  =  shear  stress  of  the  flow,  lb /ft2 
zc  =  critical  shear  stress  for  particle  movement,  lb/ft2 

The  shear  stress  of  the  flow  can  be  estimated  from  the  product  yhS,  where 
yis  the  specific  weight  of  water  (62.4  lb/fts),  h  is  the  depth  of  flow,  and  S  is 
the  surface  slope.  The  difficulty  with  the  DuBoy’s  equation,  and  other 
similar  methods  involving  the  critical  shear  stress,  is  estimating  the  critical 
shear  stress.  This  estimation  is  difficult  enough  for  pure  soil  or  sediment 
particles,  but  the  complexity  is  compounded  with  a  mixture  of  soil  and  MC 
particles,  where  the  MC  particles  can  have  much  different  particle  sizes 
and  specific  gravity.  For  example,  the  specific  gravity  of  lead  is  11.35 
compared  to  2.65  for  most  soils.  Typical  particle  sizes  for  lead  fragments 
from  a  firing  range  are  on  the  order  of  hundreds  of  micrometers  (microns, 
to-6  m)  (Larson  et  al.  2005)  compared  to  64  microns  or  less  for  fine  sand 
and  other  fine  soil  material. 

There  are  formulae  for  estimating  the  critical  flow  velocity  uc,  such  as  the 
following  equation  suggested  by  Julien  (1995)  with  uc  in  units  of  ft/sec, 

u,=lAj(G-i)gd, 


where, 

G  =  specific  gravity  of  the  particle 
g  =  acceleration  due  to  gravity,  32.2  ft2/sec 

With  the  critical  flow  velocity,  the  critical  shear  stress  can  be  estimated. 
The  approach  that  was  adopted  for  estimating  uc  for  a  mixture  was  to 
calculate  the  specific  gravity  and  average  particle  size  for  the  mixture 
taking  into  account  the  average  particle  sizes  and  the  specific  gravities  of 
the  soil  and  MC  particles  and  the  fraction  by  weight  of  the  MC  in  the  soil. 

The  Duboy’s  equation  was  tested  and  results  compared  with  those 
computed  with  the  Modified  Universal  Soil  Loss  Equation  (MUSLE)  for  soil 
only  (no  MC).  The  comparison  was  relatively  close.  Next  the  Duboy’s 
equation  was  applied  to  laboratory  studies  reported  by  Larson  et  al.  (2004). 
The  laboratory  studies  consisted  of  boxes  tilted  at  0.0625  slope  filled  with 
metal-contaminated  soil  and  exposed  to  controlled  rainfall  for  a  period  of 
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16  weeks.  Runoff  and  leachate  from  each  box  (i.e.,  lysimeter  or  test  cell) 
were  periodically  collected  and  analyzed  for  dissolved  and  total  metals,  as 
well  as  total  suspended  solids  (TSS).  The  average  total  runoff  of  soil  mass 
from  all  of  the  eight  test  cells  at  the  end  of  the  16  weeks  was  165  g.  Duboy’s 
equation  was  applied  to  the  conditions  of  the  laboratory  experiments.  The 
shear  stress  computed  for  the  laboratory  conditions  was  1.0E-4  lb/ft2,  which 
is  much  less  than  the  computed  critical  shear  stress  for  the  soil  of 
6.88  lb/ft2.  Thus,  zero  erosion  was  computed.  This  was  the  case  even  for 
pure  soil  with  no  metal  content.  Further  examination  of  the  Duboy’s 
equation  revealed  that  results  are  extremely  sensitive  to  the  value  of  the 
critical  shear  stress,  where  the  value  had  to  be  very  close  to  but  slightly  less 
than  the  value  of  the  flow  shear  stress  in  order  to  get  the  observed  amount 
of  soil  erosion.  It  was  concluded  that  a  method  was  needed  that  did  not 
require  estimation  of  the  critical  shear  stress. 

One  of  the  interesting  points  of  the  measured  laboratory  data  is  that  the 
concentration  of  lead  in  the  runoff  was  about  1.1  percent  of  the  sediment 
concentration,  which  is  close  to  the  initial  soil  concentration  of  lead,  which 
was  8100  mg/kg,  or  0.81  percent  lead.  This  could  be  a  coincidence,  but  it 
could  mean  that  lead  particles  are  eroding  as  a  mixture  with  the  sediment. 
It  could  also  mean  that  small  lead  particles  are  dissolving  and  adsorbing  to 
small  sediment  particles  that  are  being  eroded. 

Einstein  and  Brown  Equations 

The  Einstein  and  Brown  equations  do  not  require  estimation  of  a  critical 
shear  stress.  The  equations  are  not  repeated  fully  here  and  can  be  found  in 
the  text  by  Julien  (1995).  There  is  a  non-linear  relationship  for 
dimensionless  volumetric  unit  sediment  discharge  ( qbv *)  versus  the 
dimensionless  Shields  parameter  r* .  The  variable  qbv*  is  defined  as 


Qbv * 


Qbv 


(A3) 


where  co0  is  Rubey’s  clear-water  sediment  fall  velocity,  an  equation  for 
which  is  presented  by  Julien  (1995)  that  is  dependent  on  soil  grain  size 
and  specific  gravity.  The  dimensionless  Shields  parameter  is  defined  as 


T 


(ys-vK 


(A4) 
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where  ys  is  the  specific  weight  of  sediment,  and  all  other  variables  have 

been  previously  defined.  There  are  three  non-linear  equations  used  to 
relate  qbv*  to  r»  depending  on  the  value  of  r, . 

Application  of  the  Einstein  and  Brown  equations  to  the  laboratory 
conditions  yields  a  cumulative  soil  erosion  of  0.2  g,  which  is  three  orders 
of  magnitude  smaller  than  the  measured  soil  erosion  amount  of  165  g. 
However,  the  value  of  the  Shields  parameter  for  the  laboratory  study  was 
less  than  reported  field  values,  thus  outside  the  range  of  values  used  to 
develop  the  relationships. 

The  Einstein  and  Brown  equations  were  also  applied  to  conditions  at  Fort 
A.P.  Hill,  VA.  Fort  A.P.  Hill  was  a  proof-of-concept  test  case  for  the 
TREECS  Tier  1  models  (Dortch  et  al.  2010).  The  average  annual  runoff 
estimated  for  Fort  A.P.  Hill  of  0.067  m/yr  was  used  for  the  calculations  to 
test  the  Einstein  and  Brown  equations.  The  soil  erosion  computed  with  the 
Einstein  and  Brown  equations  was  30  tons/acre-yr,  compared  with 
54  tons/acre-yr  as  computed  with  the  Universal  Soil  Loss  Equation 
(USLE).  Thus,  it  appears  that  the  Einstein  and  Brown  equations  hold  some 
promise  for  computing  soil  erosion  at  field  scales. 

The  good  traits  about  the  Einstein  and  Brown  equations  are  that  they  do 
not  require  a  critical  shear  stress.  They  are  dimensionless  —  results  will 
vary  with  runoff  flow  rate  and  soil  conditions  —  and  they  include  variables 
that  could  be  used  to  describe  characteristics  of  sediment  and  MC 
mixtures.  The  particle  specific  gravity,  specific  weight,  and  size,  can  be 
adjusted  to  reflect  a  mixture.  The  fraction  by  weight  of  soil  and  MC  can  be 
used  to  adjust  these  variables  by  assuming  spherical  particles. 

Recommendations 

It  is  recommended  that  detailed  bed  load  transport  equations  for 
sediment-MC  mixtures  not  be  used  in  Tier  2  at  this  time.  If  the  user  wants 
to  consider  erosion  of  solid  phase  particles,  then  such  an  option  should  be 
included  in  Tier  2,  but  the  calculation  should  be  based  solely  on  results  of 
the  USLE  that  is  already  needed  for  soil  erosion  that  exports  aqueous 
phase  (adsorbed  and  dissolved)  MC.  This  will  mean  that  solid  phase  MC 
will  not  affect  the  overall  soil  erosion;  rather,  it  will  be  considered  as  part 
of  the  soil  erosion.  The  laboratory  results  discussed  above  indicate  that 
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this  may  not  be  a  bad  assumption.  It  is  anticipated  that  inclusion  of  this 
term  will  be  for  sensitivity  assessment  in  most  cases. 

In  this  case,  the  solid  phase  erosion  export  flux,  Fes  (g/yr),  can  be 
computed  from 


Fes=fMCPbAE 


(A5) 


where, 

f mc  =  fraction  by  weight  of  solid  phase  MC  mass  to  soil  mass 
pb  =  soil  dry  bulk  density,  g/m3 
A  =  AO  I  site  area,  m2 

E  =  soil  erosion  rate  as  determined  from  the  USLE,  m/yr 


It  is  recognized  that  /mc  can  be  computed  from 


J 'mc 


AT, 


AZ,,Ph 


(A6) 


where  Ms  (g)  is  the  solid  phase  MC  mass  in  the  AOI,  and  Zb  (m)  is  the 
thickness  of  the  surface  soil  layer.  Plugging  Equation  A6  into  A5  yields 


Equation  Ay  resembles  Equation  23  for  Fe  after  recognizing  that 
concentration  is  mass  divided  by  soil  volume  (AZb).  If  the  user  elects  to 
include  solid  phase  MC  erosion,  then  the  additional  flux  term  Fes  will  be 
computed  from  Equation  Ay  and  subtracted  from  the  right  side  of  the  solid 
phase  mass  balance  equation  as  noted  in  Equation  2.  The  export  flux  Fes 
will  be  added  to  the  particulate  mass  flux  of  the  WFF  to  import  into  the 
receiving  water  model.  If  this  option  is  not  selected,  then  Fes  will  be  set  to 
zero. 

Preliminary  testing  of  the  Tier  2  soil  model  for  a  constant  source  loading  of 
RDX  was  performed  with  Fes  included  using  Equation  Ay  to  compute  it. 
Testing  with  1,000  micron  particles  indicated  that  the  mass  export  due  to 
this  term  can  be  about  the  same  order  of  magnitude  as  the  export  due  to 
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rain  induced  pore  water  extraction  (fv)  and  roughly  two  orders  of 
magnitude  greater  than  soil  erosion  of  the  aqueous  dissolved  and  adsorbed 
phases  (Fc).  For  100  micron  particles,  mass  export  due  to  Fes  was  about 
6  times  greater  than  that  due  to  Fe,  but  export  due  to  Fr  was  about  20  times 
greater  than  that  due  to  Fes.  Thus,  including  the  Fes  term  can  be  relatively 
important  compared  with  Fe  but  less  important  compared  with  Fr.  However, 
the  relative  importance  of  the  term  depends  on  various  factors  including 
particle  size,  the  erosion  rate  E  of  the  site,  soil  partitioning  Kd,  and  water 
solubility.  As  E  increases,  the  relative  importance  of  Fes  increases. 

The  flux  Fes  is  an  additional  mass  loading  term  in  the  WFF  file  for  the 
surface  water  models  as  noted  in  Table  4.  However,  the  fact  that  this  loading 
is  solid  phase  MC  mass  will  not  be  recognized  by  the  surface  water  models 
since  these  models  do  not  handle  the  non-aqueous  phase.  This  extra  loading 
will  be  added  to  the  other  two  loadings  (Fr  and  Fe)  and  treated  as  part  of  the 
total  aqueous  phase  mass  loading.  This  means  that  solid  phase  mass  will  be 
immediately  available  for  partitioning  between  sediment  and  water  just  as 
the  non-solid  phase.  For  HE,  this  is  not  a  bad  assumption  since  dissolution 
rates  are  relatively  fast.  It  is  also  not  a  bad  assumption  for  metals,  since  they 
have  high  sediment  partitioning  coefficients  which  tie  up  aqueous  phase 
adsorbed  metals  in  the  sediments  where  the  solid  phase  metals  will  be 
located,  too.  Thus,  the  sediment  total  concentration  should  not  be  much 
different  whether  or  not  the  solid  phase  mass  is  treated  separately. 

The  TREECS  advisory  panel  expressed  an  interest  in  the  development  of 
an  additional  tool  for  assessing  single  storm  events  and  first  flush  MC 
concentrations  in  streams.  It  is  recommended  that  the  Einstein  and  Brown 
equations  be  considered  for  erosion  of  soil-MC  mixtures  in  the  transient, 
event-oriented  TREECS  soil  model  if  storm  event  capabilities  are  added  to 
TREECS. 
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Appendix  B:  Dissolution  Model  Comparisons 

Comparison  with  Outdoor  Results 

The  dissolution  formulation  (Equation  40)  for  the  TREECS  Tier  2  soil 
model  was  applied  to  experimental  results  reported  by  Taylor  et  al. 

(2009b)  for  TNT  to  gain  validation  confidence.  The  experiments  were 
conducted  outdoors  with  individual  cm-sized  chunks  of  TNT  and  Tritonal 
(an  80:20  mixture  of  TNT  and  aluminum  flakes)  exposed  to  natural 
precipitation  for  one  year.  The  11  TNT  chunks  varied  from  0.361  g  to 
1.975  g  with  an  average  mass  of  0.922  g,  and  the  five  Tritonal  chunks 
varied  from  2.162  to  5.32  g  with  an  average  mass  of  3.154  g.  Each  chunk 
rested  on  glass  frits  in  individual  glass  funnels.  Each  funnel  was  connected 
to  a  1  L  glass  bottle  beneath  the  funnel  so  that  all  precipitation  interacting 
with  the  chunk  could  be  collected  and  analyzed.  The  glass  bottles  were  set 
inside  an  insulated  wooden  box  beneath  the  funnel  to  keep  the  collection 
in  the  dark  to  prevent  photo-transformation. 

The  TREECS  dissolution  formulation  was  compared  against  the  average  of 
the  11  TNT  and  5  Tritonal  individual  chunk  dissolution  results.  The  average 
of  the  measured  dissolved  TNT  mass  was  0.021  g  for  pure  TNT  and  0.024  g 
for  Tritonal.  However,  dissolved  TNT  accounted  for  only  about  one-third  of 
the  total  mass  losses  of  TNT  from  the  chunks.  There  was  an  unaccounted 
mass  loss  of  0.037  and  0.062  g  of  TNT  for  pure  TNT  and  Tritonal, 
respectively.  The  unexpected  and  unaccounted  mass  loss  was  attributed  to 
the  formation  and  dissolution  of  red  TNT  photo-transformation  products  on 
the  surface  of  the  chunks.  Thus,  there  is  still  much  that  is  not  well-known 
regarding  TNT  fate  processes. 

Taylor  et  al.  (2009b)  propose  a  linear  drop  impingement  dissolution 
model  for  HE  as  follows, 

where, 

Mdis(t)  =  particle  mass  loss  over  time  due  to  dissolution,  g 
AF  =  particle  surface  area  factor,  dimensionless 
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Cs  =  precipitation  temperature-weighted,  average  constituent 
solubility  limit  in  water,  g/cm3 
H(t)  =  cumulative  rainfall  over  time,  cm 
do  =  initial  particle  radius,  cm 
Do,,,  =  rainfall  mass-weighted,  mean  raindrop  diameter,  cm 

If  Equation  Bi  is  applied  for  one  year,  then  H(t)  is  Pt.  Total  rainfall, 
including  equivalent  snowmelt,  for  the  year  was  122.7  cm.  Precipitation- 
weighted  average  temperature  for  the  year  was  10.5  °C,  which  resulted  in  an 
average  solubility  of  7.1E-5  g/cm3.  The  value  of  Dom  was  estimated  to  be 
0.17  cm.  The  average  of  the  linear  model  best-fit  area  factors  AF  was  2.0  and 
1.22  for  the  pure  TNT  and  Tritonal  tests,  respectively.  Using  the  average 
initial  masses  for  TNT  and  Tritonal  of  0.922  and  3.154  g,  assuming 
spherical  particles,  and  using  the  respective  solid  mass  densities  of  1.65  and 
1.87  g/cm3,  the  initial  particle  radii  of  TNT  and  Tritonal  were  computed  to 
be  0.51  and  0.74  cm,  respectively.  Using  the  above  values  in  Equation  Bi, 
the  computed  dissolution  mass  loss  of  TNT  ( Mdis )  is  0.025  and  0.028  g, 
respectively,  for  pure  TNT  and  Tritonal.  These  values  compare  fairly  closely 
with  the  measured  dissolved  TNT  mass  of  0.021  and  0.024  g  for  pure  TNT 
and  Tritonal.  It  is  noted  that  the  same  solubility  limit  was  used  for  pure  TNT 
and  Tritonal,  which  was  inferred  as  the  recommended  approach  by  Taylor 
et  al.  (2009b). 

The  dissolution  formulation  proposed  for  the  TREECS  Tier  2  soil  model  is 
stated  as 

Fdis  =  PtaMsCs  (g2) 

where, 

Fdis  =  dissolution  mass  flux,  g/yr 
Pt  =  average  annual  total  precipitation,  m/yr  or  cm/yr 
a  =  average  specific  surface  area  of  the  solid  phase  mass,  m2/g  or 
cm2/g 

Ms  =  total  solid  phase  mass  remaining,  g 
Cs  =  annual  average  of  solubility  limit  in  water,  g/m3  or  g/cm3 

Consistency  of  units  must  be  maintained  in  Equation  B2  as  noted  in  the 
definitions  of  terms. 
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Equations  Bi  and  B2  are  quite  similar  after  recognizing  that  Fdis  and  Mdis  are 
equivalent,  as  well  as  Pt  and  H,  on  an  annual  basis.  If  Cs  is  determined  from 
precipitation-weighting  rather  than  time-averaging  of  air  temperature,  then 
it  is  equivalent  to  Cs .  The  product  aMs  represents  an  average  solid  phase 
mass  surface  area.  Thus,  the  only  additional  parameters  that  are  in 
Equation  Bi  that  are  not  in  Equation  B2  are  the  fitting  parameter  AF  and 
the  raindrop  diameter.  However,  if  the  value  of  2  is  used  for  AF  as  was  done 
for  pure  TNT,  then  the  product  of  n  AF  times  the  quantity  in  parentheses  in 
Equation  Bi  is  approximately  equal  to  the  solid  phase  mass  surface  area  for 
pure  TNT.  In  this  case,  Equations  Bi  and  B2  are  really  almost  the  same. 
Equation  B2  has  some  added  appeal  of  not  requiring  a  fitting  parameter  AF 
or  raindrop  diameter.  The  raindrop  diameter  seems  like  an  unnecessary 
parameter  when  considering  that  rainfall  often  pools  on  the  ground  causing 
rain  water  to  totally  engulf  or  submerge  each  HE  particle.  Thus,  in  real 
world  settings,  the  dissolution  rate  may  be  a  consequence  of  how  rapidly 
rain  water  on  the  ground  is  replaced,  i.e.,  the  rainfall  or  precipitation  rate, 
rather  than  the  drop  diameter  and  drop  interval. 

Equation  B2  (i.e.,  Equation  40)  was  applied  to  the  average  of  the  pure  TNT 
and  Tritonal  test  results.  In  the  application  of  Equation  B2,  Ms  was 
updated  (decreased)  by  an  amount  equal  to  the  dissolution  of  TNT.  This 
was  accomplished  by  solving  Equation  2  with  the  loading  and  solid  phase 
particle  erosion  fluxes  set  to  zero.  However,  this  update  did  not  include  the 
additional  dissolution  and  mass  loss  associated  with  supposed  photo¬ 
transformation  products.  Equation  43  was  used  to  update  the  average 
particle  diameter,  and  Equation  42  was  used  to  compute  a  given  the 
particle  diameter  and  the  solid  mass  density. 

The  application  for  pure  TNT  was  rather  straightforward,  and  the  inputs 
were  the  same  as  those  cited  for  the  linear  drop  impingement  model.  The 
values  used  for  Cs  and  psm  were  71  g/m3  and  1.65  g/cm3  (1.65E6  g/m3), 
respectively.  Assuming  spherical  particles,  this  resulted  in  an  initial  particle 
diameter  of  1.022  cm  (0.01022  m)  for  an  initial  TNT  mass  of  0.922  g,  which 
resulted  in  an  initial  value  for  a  of  3.56E-4  m2/g  (3.56  cm2/g).  Equations 
B2,  2,  42,  and  43  were  coded  into  a  spreadsheet  and  solved  using  a  first- 
order,  Euler  time  integration  with  a  time  step  of  0.01  yr.  The  computed 
cumulative  TNT  dissolved  mass  after  one  year  was  0.028  g,  which  compares 
fairly  well  with  the  linear  drop  impingement  model  results  of  0.025  g  and 
the  observed  results  of  0.021  g. 
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A  model  run  was  made  where  the  remaining  mass  Ms  was  decreased  each 
time  step  by  3  times  the  dissolution  flux  computed  for  the  time  step  in  an 
attempt  to  better  represent  the  mass  loss  due  to  dissolution  of  the  supposed 
photo-transformation  products.  This  correction  decreased  the  cumulative 
TNT  dissolved  mass  after  one  year  by  only  0.001  g.  Reducing  the  time  step 
to  0.001  yr  to  improve  numerical  accuracy  decreased  the  cumulative  TNT 
dissolved  mass  after  one  year  by  only  0.0003  g-  Overall,  the  computed 
result  for  pure  TNT  compares  very  well  with  the  linear  drop  impingement 
model  and  is  reasonably  close  to  the  measured  TNT  dissolution  mass. 

The  application  to  Tritonal  was  less  straightforward  than  for  TNT.  The  input 
water  solubility  limit  is  less  certain  since  Tritonal  has  aluminum  mixed  with 
TNT.  The  average  initial  particle  diameter  was  set  to  1.48  cm  (0.0148  m) 
based  on  an  initial  mass  of  3.154  g,  spherical  particles,  and  solid  mass 
density  of  1.87  g/cm3.  The  mass  density  was  based  on  a  ratio  of  80:20  TNT 
and  aluminum  with  respective  mass  densities  of  1.65  and  2.73  g/cm3.  With 
the  initial  particle  diameter,  the  initial  a  was  computed  to  be  2.46E-4  m2/g 
(2.46  cm2/g).  Using  the  solubility  limit  of  TNT  of  71  g/nU  and  without 
correcting  the  additional  dissolution  of  photo-transformation  products,  the 
computed  cumulative  TNT  dissolved  mass  (corrected  for  80  percent  content 
of  TNT  in  Tritonal)  after  one  year  was  0.054  g ,  or  about  double  the  linear 
drop  impingement  model  and  more  than  double  the  measured  value. 

Another  model  run  for  Tritonal  was  made  assuming  that  the  total  mass 
dissolved  each  time  step  was  3  times  the  TNT  mass  dissolved  to  account 
for  the  mass  loss  due  to  dissolution  of  supposed  photo-transformation 
products.  This  correction  only  decreased  the  cumulative  TNT  dissolved 
mass  by  0.001  to  0.053  g-  Another  model  run  was  made  assuming  that  the 
solubility  limit  of  Tritonal  was  less  than  TNT  and  in  proportion  to  the 
80:20  mixture  of  TNT  and  aluminum  with  a  zero  solubility  for  aluminum. 
This  assumption  resulted  in  a  solubility  of  56.8  g/m3.  The  cumulative  TNT 
mass  dissolved  after  one  year  for  the  lower  solubility  and  the  corrected 
total  mass  dissolution  was  0.043  g,  which  is  still  nearly  double  the 
measured  value  of  0.024  g.  It  is  noted  again  that  the  parameter  AF in  the 
linear  drop  impingement  model  had  to  be  reduced  by  nearly  half  from  2  to 
1.22  for  modeling  Tritonal.  There  is  no  fitting  parameter  for  the  TREECS 
dissolution  model  other  than  trying  to  determine  an  appropriate  solubility 
for  a  mixture.  The  relatively  large  initial  mass  of  the  explosive  and  the 
solubility  limit  are  the  primary  variables  responsible  for  the  over¬ 
prediction. 
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The  approach  described  above  for  Tritonal  was  to  model  the  dissolution 
and  total  mass  of  the  mixture  using  solubility,  initial  mass,  and  solid  mass 
density  of  the  mixture.  The  problem  lies  in  knowing  the  true  controlling 
solubility  limit  of  a  mixture.  Taylor  et  al  (2009a)  found  that  for  mixtures 
(such  as  Comp  B  and  Octol),  the  constituents  of  lower  solubility  (such  as 
RDX  or  HMX)  can  control  or  lower  the  dissolution  rate  for  the  more 
soluble  constituents  (such  as  TNT).  They  did  not  observe  an  aluminum 
controlled  dissolution  of  TNT  for  Tritonal;  rather,  TNT  dissolved  at  rate 
that  was  independent  of  aluminum  in  the  mixture.  The  small  particle  drop 
impingement  model  (Taylor  et  al.  2009a)  does  not  require  the  AF 
parameter,  and  the  model  matched  measured  dissolution  mass  loss  very 
well  for  Tritonal  over  about  75  days  using  the  solubility  of  TNT. 

If  the  solubility  of  Tritonal  is  set  to  half  of  that  of  TNT,  or  35  g/m3,  the 
cumulative  TNT  mass  dissolved  after  one  year  computed  with  the  TREECS 
dissolution  model  was  0.027  g,  which  is  close  to  the  same  result  computed 
with  the  fitted  linear  drop  impingement  model.  Perhaps  the  dissolution  of 
TNT  is  eventually  impeded  by  the  presence  of  the  aluminum  flakes  after  a 
longer  rainfall  exposure  time  of  a  year,  similar  to  the  way  less  soluble  RDX 
and  HMX  impede  the  dissolution  of  TNT  in  Comp  B  and  Octol  (Taylor  et 
al.  2009a).  Regardless  of  the  reason  for  the  over-prediction  of  TNT 
dissolved  mass  for  Tritonal,  it  is  evident  that  more  study  is  required  for 
modeling  the  dissolution  of  mixtures. 

Testing  of  the  TREECS  dissolution  model  for  HE  indicated  that  although 
the  particle  specific  surface  area  a  increases  over  time  due  to  particles 
shrinking  as  dissolution  occurs,  the  actual  surface  area,  which  is  the 
product  a  Ms,  decreases  over  time  because  Ms  is  decreasing  at  a  greater 
rate  than  a  is  increasing.  Thus,  the  dissolution  rate  flux  decreases  over 
time.  This  result  is  for  the  case  on  no  additional  HE  loading. 

Comparison  with  Laboratory  Results  for  Comp  B 

The  results  of  a  laboratory  dissolution  study  of  Comp  B  by  Lever  et  al. 
(2005)  were  also  used  to  evaluate  the  proposed  TREECS  dissolution 
model.  Water  was  dropped  on  individual  particles  of  Comp  B,  and  the 
dissolved  masses  were  measured  in  the  laboratory  experiments.  The 
TREECS  model  (Equations  B2  or  40,  2,  42,  and  43)  was  applied  for  the 
conditions  of  particle  1  in  the  paper,  which  had  an  estimated  initial  Comp 
B  mass  of  1.821  mg.  The  initial  particle  diameter  was  computed  assuming 
a  sphere  and  a  particle  density  of  1.65  g/cm3. 
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With  a  rainfall  rate  of  0.55  cm/hr,  1.73  mg  of  Comp  B  was  measured  as 
dissolved  in  68  days  with  1.09  mg  of  RDX  and  0.64  mg  of  TNT  dissolved  in 
the  68  days.  Using  the  mass  ratio  of  RDX/TNT  in  the  Comp  B  samples  of 
about  1.72,  the  computed  solubility  limit  of  the  mixture  was  76.9  mg/L. 
Using  the  mixture  solubility,  the  TREECS  dissolution  model  predicted  that 
1.74  mg  of  Comp  B  would  be  dissolved  in  68  days  with  1.1  mg  of  RDX  and 
0.64  mg  of  TNT  dissolved.  The  RDX  and  TNT  dissolution  masses  are 
simply  the  Comp  B  dissolution  mass  times  the  respective  mass  fractions  of 
each  component. 

If  the  solubility  of  RDX  is  used  as  suggested  by  Lever  et  al.  (2005)  and 
Taylor  et  al.  (2009a)  to  control  the  rate  of  Comp  B  dissolution,  then  the 
TREECS  model  under-predicts  the  dissolved  total  mass  after  68  days  with 
only  1.41  mg  dissolved,  rather  than  1.74  mg.  If  the  mass  and  solubility  of 
the  individual  constituents  are  used,  then  the  dissolved  mass  of  RDX  is 
0.96  mg  after  68  days  or  is  under-predicted,  and  the  mass  of  TNT  is 
predicted  to  be  gone  in  44  days,  or  the  dissolved  mass  is  over-predicted. 
Therefore,  these  results  indicate  that  it  is  better  to  use  the  solubility  of  the 
mixture.  However,  as  discussed  in  Chapter  8,  there  are  other  factors  to 
consider  before  recommending  modeling  dissolution  of  mixtures  in 
TREECS  Tier  2. 


Evaluation  of  Dissolution  Time 


Equation  B2  was  used  to  estimate  the  length  of  time  it  would  take  to  nearly 
fully  dissolve  a  l-cm  particle  of  Comp  B.  Equation  B2  can  be  rewritten  as 


F  -M. 

1  dis  ~ 


dt 


yms 


(B3) 


where  y  =  PtaCs  and  has  units  of  yr1.  Equation  B3  can  be  solved  for  the 
time  required  to  reached  90  percent  mass  loss  via  dissolution,  resulting  in 


t 


90 


2.3 

y 


(B4) 


where  t90  (yr)  is  the  time  required  to  reached  90  percent  mass  loss. 

As  an  example,  it  is  assumed  that  a  is  constant  over  time  for  a  l-cm  chunk 
of  Comp  B  with  a  density  of  1.65  g/cm3  and  a  solubility  of  40  mg/L 
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(corresponding  to  RDX),  which  is  suggested  by  Taylor  et  al.  (2009a).  The 
chunk  is  assumed  to  be  exposed  to  rainfall  of  75  cm/yr,  which  is  the 
average  for  the  U.S.  as  suggested  by  Taylor  et  al.  (2004).  With  these 
conditions,  Equation  B4  predicts  that  90  percent  of  the  Comp  B  chunk  will 
dissolve  in  210  years.  The  time  should  be  less  than  210  years  since  a  (thus 
y)  will  increase  over  time.  In  fact,  the  full  numerical  model  (i.e.,  solution  of 
Equations  B2,  2,  42,  and  43)  predicts  that  90  percent  of  the  Comp  B  mass 
is  dissolved  in  147  years.  These  predictions  for  nearly  full  dissolution  time 
are  far  less  than  those  suggested  by  Taylor  et  al.  (2004),  which  were  on  the 
order  of  1,000  years  or  more.  Given  that  the  larger  chunks  of  HE  will  tend 
to  fragment  into  smaller  particles  over  time,  the  complete  dissolution  time 
could  be  even  far  less  than  these  projections.  As  an  example,  the  TREECS 
model  predicts  that  a  1  mm  particle  of  Comp  B  exposed  to  75  cm/yr  of 
rainfall  will  be  dissolved  90  percent  in  15  years.  Compared  to  less  soluble 
MC,  such  as  metals,  HE  particles  dissolve  rapidly. 

Taylor  et  al.  (2009b)  used  their  more  recent  linear  drop  impingement 
model  to  project  long-term  dissolution  and  life  spans  of  various  sized  HE 
particles.  Their  results  show  that  it  would  take  about  too  years  to  fully 
dissolve  1  g  of  TNT  at  an  annual  average  temperature  of  10.5  °C  with  an 
average  annual  rainfall  of  too  cm/yr.  Equation  B4  predicts  it  would  take 
94  years  to  dissolve  90  percent  of  the  1  g  chunk  of  TNT.  The  TREECS 
numerical  dissolution  model  predicts  it  will  take  65  years  to  dissolve 
90  percent  of  the  1  g  chunk  of  TNT  and  95  years  to  dissolve  99  percent  of 
it.  This  latter  result  compares  closely  with  the  too  year  projection  made  by 
Taylor  et  al.  (2009b).  As  they  pointed  out  and  as  can  be  observed  from 
Equations  B2  and  B4,  dissolution  rate  and  particle  life  span  are  linearly 
related  to  precipitation  rate.  Thus,  if  the  precipitation  is  half  of  the  above 
rate,  the  life  span  will  double. 
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